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TABLE OF SYMBOLS 



Symbol 




Quantity 


Units 


A 

m 


= 


cross sectional area of metal 


ft 2 


A 

f 


= 


cross section area of fluid 


ft 2 


a 


= 


inner radius of pipe wall 


ft 


a 


= 


half width of channel (slab geometry) 


ft 


a 


= 


^1^2^ * usec * equations (IV-29) 

and (IV-31) 


hr" 1 


B f 


= 


^ a /Kf = Fluid Biot number defined 








by equation (V-33) 




Bf* 


= 


2ha./K^ _ Fluid Biot number defined 
by equation (V-62) 




B 

m 


= 


kb/K = pipe wall Biot number 
defined by equation (V-33) 




b 


= 


outer radius of pipe 


ft 


b 


= 


thickness of wall (slab geometry) 


ft 


b i 


= 


/ j^CjA^ = parameter defined by 
equation (IV- 2) 


hr * 


b 2 


= 


hP 

/ / > m C m A m = parameter defined 
by equation (IV_5) 


hr' 1 


°f 


- 


Fluid specific heat 


BtU /lbm°F 


c m 


= 


pipe wall specific heat 




Ci 




coefficient defined by equation (IV-57), 
i = 0, 1 




c * 

i 


= 


coefficient defined by equation (IV-67), 
i = 0, 1 
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h 

h 

h* 

h.((, n;t) 
X n 

J n (x) 

K ma 

Kfflt 

K 

m 

Kfa 

K ft 

K 

f 

k i 

V 

K n (x) 

P 

q 

r 

s 



= heat transfer coefficient at fluid -wall 
interface 

= heat transfer coefficient defined by- 
equation (II- 15) 

= heat transfer coefficient defined by- 
equation (II- 9) 

= Function defined by equation (IV- 45) 

= modified Bessel function of first kind 
of order n 

= complex integral defined in Appendix A 
equation (5) 

= Bessel function of first kind of order n 

= pipe wall axial conductivity 

= pipe wall transverse conductivity 

= pipe wall isotropic conductivity 

= fluid axial conductivity 

= fluid transverse conductivity 

= fluid isotropic conductivity 

= - coefficients defined by equation 
(IV-40), i = 0, 1,2, 3, 4 



Btu /ft 2 hr°F 

Btu /ft 2 hr°F 

Btu/ft 2 hr°F 



Btu/ft hr°F 
Btu/ft. hr°F 
Btu/ft. hr°F 
Btu/ft. hr°F 
Btu/hr ft °F 

Btu /hr ft °F 



= coefficients defined by equation 

(IV-42), i = 1,2, 3, 4, 5 

= modified Bessel function of second 
kind of order n. 

= perimeter 

= (s/ * ) 

= radial coordinate 
= Laplace transform parameter 



ft 



ft 



hr 



"I 
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T 

T* 

t 

u. 

1 

V 

W(x,y,t) 
W(r, x, t) 
x 

y 

Y n (x) 



c 

1 

e 

6 * 

I 

« 

r 

t 



pipe wall temperature 

radially averaged pipe wall 
temperature 

Time 

temperature functions used in 
Section III 

fluid velocity 

temperature function defined by 
equation ( III- 21) 

temperature function defined by 
equation ( 111 - 92 ) 

axial coordinate 

transverse coordinate 

Bessel function of second kind of 
order n 



°F 

°F 

hr 

°F 

ft/hr 

°F 

°F 

ft 

ft 



GREEK LETTERS 



fluid thermal diffusivity 
pipe wall thermal diffusivity 
V (used in Section IV) 
x/V (used in Section IV) 
fluid temperature 

radially averaged fluid temperature 
( 1 ~t)/ 4 (used in Section IV) 
fluid density 
H_-t (used in Section IV) 
fluid temperature function at interface 



ft 2 /hr 

ft 2 /hr 

hr /ft 

hr 

°F 

°F 

ft 

lbrn/ft^ 

hr 
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X 



function defined by equation (V-43) 



Y 



pipe wall temperature function at 
interface 



°F 



Also, additional notations are occasionally introduced for 
temporary use in the analytical developments which follow. In such 
cases, the new notations will be adequately defined in the text. This 
also implies that the same symbol may have different significance in 
different portions of the text. There should be no confusion arising 
from this practice. In particular, different notation is frequently 
introduced for dummy variables of integration. Also, the symbols 



in the foregoing list. 

We also remark that equations are numbered serially in each 
section. In referring to equations, the serial number alone will 
indicate the correspondingly numbered equation in the same section 
as the reference; otherwise the section number will also appear. Thus, 
in Section V, for example, the reference (12) is to equation (V-12) 
and the reference (II-7) is to equation (II-7), that is, to equation 7 of 
Section II. A reference such as IID2 is to the text of subsection IID2 
and not to an equation. 




% n > \ n > JU f etc. , used as eigen values, do not appear 
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ABSTRACT 



Transient temperature distribution in axisymmetric systems of 
incompressible fluid flow through semi-infinite, externally insulated 
cylindrical pipes is considered under the assumption of constant uniform 
velocity and constant material properties. The general problem, includ- 
ing all significant heat transfer mechanisms, except radiation, is 
reduced to an integral equation containing only one unknown function, 
namely an interface temperature as a function of time and axial position. 
By neglecting one or more heat transfer mechanisms, several other 
formulations are obtained. These are analyzed so as to provide 
approximate and closed form solutions not previously available. 

Similar formulations and results are obtained for the case of fluid 
between parallel plates. 



9 



ACKNOWLEDGEMENT 



The author would like to take this opportunity to thank Professor 
J. E. Brock, Department of Mechanical Engineering, Naval Post- 
graduate School, for all the help and guidance that he so kindly 
provided for accomplishing this work. His enthusiasm and interest 
during the course of this project surpassed what might have been 
expected. Our coordinated work on this project has led to the develop- 
ment of a friendship beyond that of a student and advisor, which I shall 
always cherish. 

The help of Professor C. Comstock, Department of Mathematics, 
Naval Postgraduate School, is also gratefully acknowledged. He pro- 
vided constructive advice throughout and contributed much of the 
analysis relating to the singular perturbation development, Section IV. 

I would also like to thank Professor M. D. Kelleher, Department 
of Mechanical Engineering, Naval Postgraduate School, for reading 
through the manuscript and providing valuable suggestions. 

Lastly I would express my thanks to my wife and children who 
sacrificed so much during this time and provided me with the inspira- 
tion and enthusiasm that I needed. 



10 



I. INTRODUCTION 



A, PRELIMINARY REMARKS 

The problem of transient temperature fields in systems of fluid 
flow through circular pipes or through channels formed by parallel 
plates is of considerable practical and theoretical importance. On the 
practical side, one may consider, for example, the problem of the 
calculation of thermal stresses in main coolant piping systems of 
nuclear power plants in which abrupt changes in power output may 
cause very rapid variation in the temperature of the coolant fluid 
leaving the reactor. In such cases one may be concerned with the time 
history of the coolant temperature at some critical point, such as a 
pump or valve, downstream from the reactor. 

From a theoretical point of view, the problems implied by such 
pipe flows are of great interest. They involve two dependent variables, 
the temperatures of the fluid and of the pipe wall, each of which is a 
function of three independent variables: time and two space coordinates. 

These functions must satisfy appropriate partial differential equations 
and external boundary and initial conditions and must also conform to 
each other in an appropriate way at their common interface. Such 
problems are called "conjugated. " Few solutions for such problems 
are given in the literature. Those which were previously known usually 
involve drastic simplifications such as considering only steady state 
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behavior or restricting space variation to one dimension, or considering 
the pipe to have infinite thickness, or neglecting important mechanisms 
of heat transfer, or a combination of these. 

Accordingly it is important to consider such problems in as gen- 
eralized a form as possible both for the purpose of obtaining compact 
alternate formulations and also for the purpose of obtaining closed 
form and approximate solutions to problems in which the degree of 
simplification is less than that involved in previously known solutions. 

It is believed that the work reported herein is successful in both 
regards. 

B. STATEMENT OF THE GENERAL PROBLEM 

It is appropriate at this point to make an explicit statement of the 
general problem which will be considered herein. We are concerned 
with the flow of an incompressible fluid axially down the length of a 
semi-infinite cylindrical pipe whose cross section is a circular annulus 
of constant finite wall thickness. It is presumed that axisymmetrical 
conditions prevail so that the space variables are radius r and axial 
position x from the accessible end of the pipe. The outside of the 
pipe is presumed to be perfectly insulated so that there is no transfer 
of heat through this surface. (This condition is approximately satisfied 
in many cases of technical importance. ) 

It is assumed that the physical properties of both the fluid and the 
pipe wall are constant -- that is, they do not vary with respect to time, 
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temperature, or position. We presume that the thermal conductivity 
of the materials involved is orthotropic; that is, we provide for 
different values of conductivity in the radial and axial directions. This 
permits later specialization in different applications, either to the case 
of isotropic conductivity or to the cases where one or more conducting 
mechanisms are suppressed. Subject to this introduction of orthotropic 
conductivity, it is assumed that Fourier's law of heat conduction holds 
in each medium. In most of the work done herein, we provide for an 
interface surface heat transfer coefficient to account for surface 
resistance or other effects. This coefficient may be taken as infinite 
if there is no such resistance. However, in all cases it is assumed to 
be a constant which does not depend on time, position, or temperature. 

We have assumed the nature of the fluid flow, namely, we have 
assumed that the velocity of each element of fluid is the same, a con- 
stant velocity V in the axial direction x. This is frequently called 
slug flow. It is recognized that this is a particular case which is at 
variance with the flow conditions which would actually obtain under 
any circumstances. However, it is not far removed from the flow 
conditions which would prevail with fully turbulent flow and it does 
provide a necessary simplification which has permitted us to obtain 
useful results. 

We have confined attention to one case of temperature variation 
which can be described as follows. It is presumed that the fluid has 
been entering the pipe at zero (or reference) temperature for a 
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solution for the case of cylindrical geometry. In fact, in some cases, 
when certain simplifying assumptions and procedures are introduced, 
the two cases may become one and the same. 

C. SURVEY OF LITERATURE 

The study of heat transfer problems in systems of fluid flow 
through pipes dates back to 1883 when one of the most celebrated 
problems of heat transfer was solved by L. Graetz [11].* In Graetz's 
problem the velocity field is taken to be known, and heat transfer 
takes place under steady state conditions with pipe wall temperature 
distribution specified a priori. In recent years other problems of 
fluid flow have been treated by many authors. These studies include 
such effects as heat transfer by axial conduction and mechanical heat 
generation due to friction, i. e. , viscous dissipation. There is an 
abundance of literature in this area of heat transfer engineering. 
However, these problems are concerned solely with behavior of the 
fluid; the action of the containment is considered only in the specifica- 
tion of boundary conditions for the fluid. 

However, the transient case of temperature field in systems of 
flow through pipes or ducts when the wall thermal state is not specified 
a priori, and it is required that its effect on the temperature field in 
the fluid be accounted for, clearly falkin the class of conjugated 

^Numbers in square brackets refer to the Bibliography, page 154 
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problems of heat transfer. These conjugated problems are more 
involved inasmuch as the mutual thermal effects of the wall and the 
fluid must be accounted for. Accordingly their treatment presents 
serious mathematical difficulties. Even when the velocity field is 
assumed to be known, the thermal problems are still of great difficulty, 
mainly because of the coupling that is introduced by the condition at 
the common boundary of the fluid and the wall. 

At present there are numerous publications dealing with conjugated 
problems of heat transfer. In all cases, however, they seem to be 
particularized and as yet no general method has been obtained for 
treatment of these problems. 

An early treatment of conjugate problems of heat transfer was 

* 

given by T. L. Perelman [25] who considered steady state heat transfer 
to a fluid flowing around a body containing internal heat sources. He 
also considered the asymptotic solution to the type of integral equation 
that occurs in the analysis of such conjugated problems: namely, a 

Fredholm equation of second kind with singular kernel and infinite 
range. He considered slip flow around a body and also considered 
the laminar boundary layer flow over a thin plate. Perelman neglected 
the axial heat conduction in the fluid. 

N. A. Avdonin [4] considered transient behavior in a similar case 
but with no heat sources present. He also neglected the axial conduc- 
tion in the fluid. He solved the problem by using Green's functions 
and Laplace transforms. He treated the conjugated problem by 
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reducing it to a Volterra integral equation which he solved using 
Laplace transforms. 

D. E. James and N. G. William [ 1 7] considered the steady state 

problems of Poiseuille - Couette flow between parallel planes. In 

this work the boundary planes were taken to be of finite length in the 

direction of flow and of finite thickness perpendicular to the direction 

of flow. These authors assumed an arbitrary temperature distribution 

at the fluid-wall interface and solved the energy equations separately 

for each material for a constant interface temperature distribution. 

Then they used Duhamel’s superposition integral to obtain the general 

/ 

solutions in terms of the arbitrary interface temperature distribution. 
Next applying the continuity condition for the heat flux they obtained an 
integro-differential equation for the interface temperature distribution. 
They solved this by the method of collocation. 

Recently A. V. Luikov, V. A. Aleksashenko, and A. A. 
Aleksashenko [22], considered conjugated steady state heat transfer 
in circular and plane tubes for the case of developed Poiseuille velocity 
distribution allowing for mechanical energy dissipation. They state 
that general methods have been developed for both steady and unsteady 
heat transfer, but they give no example or discussion for other than 
steady state cases. They essentially reduce the problem to a singular 
integral equation for the unknown temperature at the fluid -wall 
interface. This is a very interesting paper. 
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E. M. Sparrow and F. N. DeFarias [33] investigated the laminar 
flow through parallel plate channels under the assumption of sinu- 
soidally varying fluid inlet temperature. This work is in the spirit of 
the new results to be reported herein inasmuch as a solution for 
sinusoidally varying end temperature can be related to that for step 
change in end temperature and conversely. Accordingly it is important 
to delineate clearly what simplifying assumptions they employed. They 
neglected axial conduction in both fluid and pipe and they assume pipe 
temperature varies with axial position but not with transverse position. 
They considered constant fluid velocity as we do herein. It happens 
that these simplifying assumptions do not coincide with any set 
treated in this thesis. 

V. S. Arpaci [3] considered a transient problem of the type in 
which we are interested. Assuming as we do, the case of slug flow, 
he uses the term "rod" to describe what we call the fluid. He shifts 
viewpoint in an obvious and nonessential way by considering the pipe 
to be moving at constant velocity with its contents at zero velocity. 
Axial conductivity is neglected in both substances. Radial conductivity 
is accounted for in the fluid* Radial conductivity of the pipe wall 
is taken as infinite. The transient effect is due to sudden and uniform 
heat generation in the contents. Like ours, his pipe has perfect 
external insulation. He used the method of finite Hankel transforms 
to reduce the problem to .the solution of an integro -differential equation 
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which was solved by Laplace transforms. For large values of time an 
asymptotic behavior of the temperature function was found, again using 
Laplace transforms. 

Other significant contributions to the study of conjugated heat 
transfer problems have been made by M. D. Kelleher and K. T. Yang 
[19], R. Siegel [34], and R. Siegel and J. M. Savino [35]. 

F. E. Moreland [23] employed numerical inversion of Laplace 
transformation to deal with heat exchanger problems which are similar 
to our problems, except for suppression of all heat diffusion in the fluid. 
D. D. Hiep [13] extended this work significantly by including the effect 
of axial heat diffusion. (His treatment is for porous media, but our 
formulation could be similarly generalized by use of some fictitious 
material and geometrical properties. ) His study contains a large 
bibliography. The problem considered in Section VI hereof is a special 
case of that considered by Hiep. We are able to obtain analytical 
results whereas Hiep*s analysis is purely numerical. 

Two earlier solutions which provide significant points of departure 
should also be mentioned, even though idealizing assumptions made 
the analyses very simple. One of these is by W. Munk [24] who 
assumes identical constant temperature in both pipe and contents at 
any cross section. He neglects all diffusion of heat and arrives at a 
solution which repeats the time-temperature history at the entrance 
point, but with a propagation velocity less than the velocity of the 
fluid. The second, to which many investigators have made contributions, 
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is discussed in section 15. 3 of [7] where many references are cited. 
Arpaci [2], page 353 et seq. , gives a very compact treatment of this 
problem in which there is no heat diffusion in either the axial or radial 
direction. 

D. SUMMARY AND OUTLINE OF SUCCEEDING SECTIONS 

The following telegraphically brief statements are intended to 
provide a summary of the content of the succeeding sections. 

Section II - Derives basic equations which will be employed in 
subsequent sections. 

Section III - Considers radial and axial heat diffusion in each 
material. Arrives at integral equation (III-70) for slab geometry and 
integral equation (III- 181) for cylindrical geometry. No solution of 
these equations is attempted. 

Section IV _ Radially averaged temperatures are considered in 
each material. An ordinary perturbation solution is constructed 
which is valid "behind the (temperature) front. M A singular perturba- 
tion solution is developed which gives variation ahead of the front. A 
scheme for matching these at the front is given. 

Section V - By averaging the fluid temperature in the transverse 
direction and neglecting axial conduction in the pipe wall^closed form 
solutions are obtained in the case of slab geometry: equations (V-39, 

40). The case of cylindrical geometry turns out to differ only as 

regards the definition of a certain auxiliary function so that the 
solution has the same expression as for the case of slab geometry. 
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Section VI - By taking radial averages for both fluid and solid 
temperatures and by neglecting axial conduction in the solid, a 
simplified version of the preceding problem is obtained. (This may 
also be obtained by neglecting axial conduction in the solid in the case 
given in section IV). A closed form solution is given: equation (VI-21) 

and (VI-23). 

Section VII - Here the problem is greatly simplified by assuming 
(in effect) that there is no radial variation in temperature. Axial 
conduction in both materials is considered, however; (i. e. , we take 
average fluid temperature equal to average solid temperature at any 
section and employ the equations of section IV). The resulting 
solution, (VII- 10) is not new but is included here for completeness and 
comparison purposes. 

Section VIII - This section includes discus sion and recommendation 
for further study. 
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II. FUNDAMENTAL EQUATIONS 



A. GENERAL REMARKS 

All the analysis presented in this thesis relates to a definite 
physical problem which will be described shortly. The real physical 
problem is so complicated as to preclude the possibility of dealing 
with it without making certain idealizations. In what follows immediately, 
the problem will be described in mathematical terms and the idealiza- 
tions which are common to all the subsequent analysis herein will be 
delineated. The simplified problem which results is still of great 
complexity. It will be further discussed in section III. No solution of 
this problem will be presented but it will be shown that the problem 
can be cast in a more compact manner which gives some promise of 
being amenable to analysis. 

In order to obtain additional results, it has been found necessary 
to make additional idealizations or simplifications. These lead to 
several simplified problems which are specifically delineated and 
dealt with in sections IV through VII. For all of this work, however, 
the derivation given immediately below and the equations which result 
therefrom are fundamental. 

1. Physical Model 

The basic physical problem of interest is that of heat transfer 
taking place between a conduit and the fluid flowing therein. The former 
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is perfectly insulated on its exterior surface. It may be either a 
circular pipe with uniform wall thickness or a duct composed of two 
parallel plates. (In what follows these two cases will be referred to 
as "cylindrical geometry" and "slab geometry, " respectively. It is 
assumed that the fluid entering the conduit has had constant temperature 
for sufficient time that the conduit and its fluid contents have reached 
this constant temperature throughout. Then the temperature of the 
fluid entering the conduit is suddenly changed by a finite amount and 
is thereafter maintained at the new temperature. (The end surfaces 
of the conduit also experience the same step change in temperature. ) 

It is required to predict the unsteady temperature fields in both the 
fluid and the conduit for all locations at any subsequent time. 

Figure 1 illustrates the two geometries which are considered. 

That for the cylindrical geometry is self-explanatory, except to remark 
that there is no variation whatsoever with respect to angular coordinate. 
The case of slab geometry involves identical parallel plates which are 
s emi- infinite in extent in the x direction, which is the direction of 
fluid flow, and are infinite in the z direction. It is assumed that there 
is no variation whatsoever with respect to the coordinate z of any 
quantity. In making a mathematical formulation of the problem, it 
will be convenient to take advantage of the symmetry with respect to 
the plane y = 0; for this reason the auxiliary coordinates y^ and y^ 
are used as indicated. 
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a - cylindrical geometry 



b - slab geometry 



Figure I : 




Geometry and Dimensions of Conduits 
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2. Assumptions 



a. The fluid is considered to be incompressible. 

b. The fluid is assumed to be flowing with constant velocity V 
(slug flow). 

c. Material properties are taken to be constant and independ- 
ent of time, position, and temperature. 

d. It is assumed that there may be an additional thermal 
resistance (possibly due to a very thin layer of foreign material at the 
interface) the heat transfer coefficient of which is a constant. 

e. Viscous dissipation is assumed negligible. 

f. Heat transfer by radiation is considered to be negligible. 

g. Fourier's law of heat conduction is taken to be valid 
throughout. 

h. The diffusivity appearing in Fourier's law for fluid is 
assumed to include not only the molecular diffusivity but also an eddy 
diffusivity which is effective. The combined diffusivity is assumed 
to satisfy assumption c. 

i. There is variation only with respect to time, the axial 
dimension, and a radial dimension. Figure la^or a transverse 
direction, Figure lb. 
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B. DIFFERENTIAL EQUATION 



1 . Cylindrical Geometry 

There is no difficulty in establishing, by customary methods 
of derivation, the following partial differential equations: (1) which 

describes the fluid and (2) which describes the solid pipe wall. We 
have anticipated simplifications which will be made in later sections 
hereof by considering artificial media in which thermal conductivities 
in the axial and transverse directions may differ. These conductivities 
are distinguished by a second subscript, a to denote axial or t to denote 
transverse. If we are considering isotropic materials (as we do in 
section III hereof, for example) we take the axial conductivity to be 
equal to the transverse conductivity. However, if in order to make 
mathematical headway in a certain analysis we neglect axial conductivity 
in the solid (say), we can obtain the corresponding mathematical des- 
cription simply by setting K =0. 

ma 

There should be no difficulty in recognizing the significance 
of the various terms which appear in equations (1) and (2) below. 




( 1 ) 




( 2 ) 
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One may observe that (2) can be obtained from (1) by taking 
V = 0, replacing 9 by T, and replacing subscript f by subscript m. 

2. Slab Geometry 

The partial differential equations governing the temperature 
distribution in the flowing fluid and in the solid duct wall are 



K 






(3) 



for the fluid, and 



k ir , k "ir 

ma ¥ mt 




( (4) 



for the pipe wall. Observe the use of subscripts 1 and 2 for the 
coordinate y, consistent with figure lb. 



C. INITIAL AND BOUNDARY CONDITIONS 
1. Cylindrical Geometry 
a. For Fluid 

(1) Initial Condition, t = 0. The initial temperature of 
the moving fluid is considered to be uniform, and without any loss of 
generality, it may be considered as zero. 

0 = o o * x < oo , o 4 r ^ & , t = o 

(2) Boundary Condition at x = 0. At this boundary the 
entering fluid temperature is given a sudden step change to a higher 



27 



value thereafter. Without loss of generality, the amount of increase 
may be taken to be unity. 

0 = 1 *-° 

(3) Boundary Conditions as x — »qq. As x approaches 
infinity the fluid temperature decreases to its initial value. 



0 —+o X — > oo j o<f<a , 

(usually, this condition is employed in the form: 9 remains bounded 

as x — ►oo . ) 

(4) Boundary Condition at r = 0. It is required that as 
r — 0 the fluid temperature remain finite. This condition may also 
be stated by requiring that the temperature distribution in the fluid 
to be symmetric about the center line. 

Either 






oo 



o <X <°° ■> i: % 0 



r—> o 



or 



T>e 

'dr 



= o 



r= o j o <x <°° j ^ = ° 



b. For the Pipe Wall 

(1) Initial Condition, t = 0 . The initial temperature of 
the pipe wall is considered to be uniform and equal to the fluid tem- 
perature at the initial time. 



T = o o <"x < j <z Kr < b j t - o 
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(2) Boundary Condition at x = 0. The pipe wall is 



assumed to experience the same unit step change of temperature as 
that of the fluid. 

T = 1 x = o , a^r^b , t ^ o 

(3) Boundary Condition as x — * ►oo. Similarly it is 
required that the pipe wall temperature approach its initial temperature 
as x approaches infinity. 



T — v o x — ► oo, a < r < b , t :> o 



(4) Boundary Condition at r = b. Since it is assumed 
that the pipe wall is perfectly insulated at the outside surface, then 
the Fourier heat conduction law requires that 

£T = Q o ^ x < oo , r=b , t > o 

br 

c. For the Interface Between Fluid and Wall (r “ a) 

Continuity of heat flux is required at the interface between 
fluid and pipe wall. 

Contemplating a constant heat transfer coefficient h at 
this boundary, we may thus write 




K mt 



dr 

£ 

dr 



= h ( 0 - T) 



O 4 : X < CO 



= h ( 9 - T) 

J 



r = a , o 
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It is well to note at this point that if the heat transfer coefficient h 
approaches infinity, then the above conditions reduce to equality of 
temperature and heat flux, in the following manner. 



0 <00 , fs a J t > o 



2. Slab Geometry 

The details and physical significance of the choice of initial 
and boundary conditions for this case are similar to those discussed 
for the case of cylindrical geometry. We will therefore state them 
mathematically as follows. 

a. For the Fluid 

(1 ) Initial Condition t = 0 . 

0-o o <X < oo j ° * %\ ~ a > £ >o 

(2) Boundary Condition at x = 0. 

0-i X « O J a < ^ j t> O 

(3) Boundary Condition as x — > oo . 

9 s 0 X — ► oo J o<4 <a j t>o 

(4) Boundary Condition at y^ = a. 



e = t 

< M = K I 



be 

>3, 



o 



o5x < 0O , 



3.* a > 



i > o 
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(Cf. figure lb; the symmetry about the face = a provides this 
condition. ) 

b. For the Wall 

(1) Initial Condition, t = 0. 

T=0 o < x < oo , o < ^ b, t = o 

(2) Boundary Condition at x = 0. 

T = 1 x = 0 , 0-^y2^b , t > o 

(3) Boundary Condition as x >-oo. 

T = 0 x — *- oo , 0 < y 2 b , t > o 

(4) Boundary Condition at y? - b. 



= 0 o < x < oo , y-> = b , t> o 

c. For the Interface Between Fluid and Wall y^ = = 0 



■K . AL = h (0 - T) 

mt ^ 

K f i. = h (9 - T) 



o ^ x < oo , Y \ ~ Y2 ~ ° * t = ° 



* 3 . 



Observe that the algebraic signs given here are a con- 
sequence of the fact that y^ and y^ are opposite in direction; Cf* 
Figure lb. Again it is easy to observe that as h approaches infinity 
the interface condition above may be written as 



0 = T 



> o ^ x <oo , Y l =Y z =° it^o 



-K = K — 

mt 1} X ft ^1 
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D. SPECIALIZATION OF ENERGY EQUATIONS 



In what follows particularizations of equations (1) and (2) for the 
cylindrical geometry and equations (3) and (4) for the slab geometry 
are introduced, which involve various techniques for simplifying the 
energy equations. These methods involve radially averaging the fluid 
and the conduit wall temperatures by integrating with respect to radius 
(for cylindrical geometry) and with respect to the transverse co- 
ordinates (for slab geometry). Thus appropriate definitions of average 
temperatures and surface film heat transfer coefficient are introduced. 
Further simplifying assumptions are made, such as neglecting the 
axial conduction in the conduit wall; these simplifications give rise to 
interesting problems. Also since the difference between cylindrical 
and slab geometry disappears when both the fluid and the wall tem- 
peratures are averaged, only the cylindrical geometry is considered 
in this case. 

1. Cylindrical Geometry 



a. 



Radial Averaging of Temperatures 



We consider isotropic conductivity in both media: 




(5) 




( 6 ) 



We also define the average temperatures 
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and 




c*. 






( 7 ) 



( 8 ) 



and we introduce an "equivalent 11 surface film coefficient 




e*- T* 



(9) 



Note that h* relates heat flux to difference between average 
fluid temperature and average wall temperature. 

Then, by integrating equations (1) and (2) with respect to 
radius, between appropriate limits, and introducing the definitions 
above, we find 



x * * 

v ^0* i de* P i .* -r*) 



( 10 ) 



for the fluid mean temperature, and 

Sr* =_L U*. ilf 

** 2 ** Vt / 



( 11 ) 
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for the pipe wall mean temperature. Equations (10) and (11) are 
essentially those considered by Hiep [13]; although his equations 
appear to provide for more generality (in order to apply to porous 
media), a change of variable shows that they are mathematically the 
same. The introduction of h‘ may legitimately be objected to on the 
grounds that physically it would be quite difficult to evaluate. However, 
its use provides a degree of generalization compared to the physically 
reasonable simplifying assumption that radial variation of temperature 
in the fluid is so small that 0'*' = 9 (r = a) = T (r = a). Moreover, 
even in cases where h — * oo, the use of h‘ < oo permits taking into 
account some resistance to radial heat flow in the fluid. 

The initial and boundary conditions for (10) and (11) follow 
naturally from those of (1) and (2). 

(1) Initial Condition t = 0. 

0 = T" = o o ^ x < oo , t = o 

(2) Bouhdary Condition at x = 0. 

0* = T '' = 1 x = o , t > o 

(3) Boundary Condition as x — >oo. 

0 - — ► 0 , T - — > o x — ► oo , t > o 

b. Case where h'*' Approaches Infinity 

An interesting special case may be obtained if it is assumed 
that radial diffusion is so great that the temperature is essentially 
constant throughout any particular cross section. This is accomplished 
by letting h‘ approach infinity. Then from (11) we obtain 
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(12) 



( B * T *) - A~K~ ( I DT* 

^ ~ !*p ^ ) 



Substituting for (0 - T*) from (12) above in (10), and noting that as IT' 



approaches infinity, then 0 = T' 1 ', yields 



* /, . K-A» \ -tr* v at*, f. g A ~ c ~ X aT J 



(13) 



for the common fluid and conduit wall temperature distribution. The 
initial and boundary conditions for (13) follow naturally from those of 
(10) and (11). 

(1) Initial Condition, t = 0. 

T* = 0 o < X <<*> , t = 0 

(2) Boundary Condition at x = 0. 

T* = 1 x = 0 , t > 0 

(3) Boundary Condition as x — >0 o> 

T '*' — > 0 x — Voo , t > 0 



The solution of this system is treated in section VII. 

c. Case where fC = 0 

m 

By neglecting axial conduction in the solid pipe wall, a 
tractable system is obtained. We rewrite equation (10) with no change 



-'t 12 * = -L i!? t g*-T *) 

•b X 1 "bX ^ £)t Af ^ ' 



( 10 ) 
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Multiplying equation (11) by K m and allowing 
yields 

JT*- ( e*- T*) 

7>t ~ P A c 

m m 



a PP roac h zero, 



(14) 



The initial and boundary conditions are 

(1) Initial Condition at t = 0. 

9* = T* = 0 0 < oo j t 

(2) Boundary Condition at x = 0, 

0 * = 1 ° ^ O 

(3) Boundary Condition as x — »qq . 




The solution of this system is dealt with in section VI . 
d. Another Special Case 

If radial averaging is done only for the fluid and variation 
with respect to radius is retained in the solid pipe wall an interesting 
problem results. However, this situation leads to an integro- 
differential equation the importance of which does not seem to be 
commensurate with its difficulty. Accordingly an additional simplifica- 
tion is made: axial conduction is neglected in the solid pipe wall. To 

accomplish this we set K ma = 0 in equation (2) and we also write 
Kjnt = ^m* For the fluid we write Kf a = Kfj- = Kf. Note that the 
radial effect is absorbed in the radial averaging for the fluid. It is 
necessary to define a new equivalent heat transfer coefficient £ as 
follows. 
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-K 




.-K = Me’- T) *t r-.a 

F 



(15) 



Then equations (1) and (2) with 9' still defined by (7) reduce to 



ie* v: \ d * (t) 

"ix 1 «f T>* <* f it AfK f { 




(16) 



and 




1 *T 



(17) 



(17) is exactly the form taken by (2) with K ma = 0. (16) is the same 

as 10) except for the new definition of surface film coefficient. The 
equivalent heat transfer coefficient h , defined by (15) relates inter- 
face heat flux to the difference between conduit surface temperature 
and fluid mean temperature; compare with h 1 ' defined by (9). Clearly 
h V h . The reasons for introducing h are covered in the remarks 
following (11). 

The appropriate boundary and initial conditions are: 

(1) Initial Condition. The initial conditions for equations 
(16) and (17) are exactly the same as those for equations (1) and (2) 



0 < X < <X> s t ~ O 

* 



T = 0 



a <r <t> s t = ° 

* s 
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(2) Boundary Condition at x = 0. The same conditions as 



stated for equations (1) still apply to equation (16): 

e* = i x ~ o j t >o 

We observe that since equation (17) is explicitly 
independent of x, it is therefore not required to specify any condition 
at x = 0 for equation (17). 

(3) Boundary Condition as x— »>co. It is again required 
that O' 1 ' remain finite; viz. 

6 * = 0 % ^ CO y t > o 

(4) Boundary Condition at r = b. Assuming perfect 
insulation at the outer surface of the cylinder, then 



vr 






ST O 



fr b y ^ ° 



(5) Boundary Condition at r = a. At the interface of 
fluid and pipe wall solid the continuity of heat flux requires 

-K m ^I = l (e*- T) r=a J t>o 

or 

We note here that if h approaches infinity, then 
the condition of continuity of heat flux above reduces to equality of 
temperature at the fluid and wall interface. 

This system is analyzed in VB1. See discussion at 
end of IID2 for explanation of why the system of IID2 is treated first. 
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2. Slab Geometry 



For all the simplified cases in Section IID1 in which average 
temperatures in both fluid and pipe wall are considered, the geometry 



Thus, specializing the case of slab geometry in the same ways leads, 
as one may easily verify, to exactly the same equations. However, 
the case treated in IIDld still retains radial variation in the pipe wall. 
Accordingly, we expect that the slab geometry counterpart will be 
described by slightly different equations, which we now proceed to 
obtain. 

a. Differential Equation 

Treating equations (3) and (4) for the case of flat wall 
duct similarly we define 

a 



0 

Then considering isotropic conductivity in the fluid, 
integration of (3) yields 



of the situation is fully described by the parameters A^, A^, and P. 





(19) 



for the fluid. Similarly with the assumption of zero axial conductivity 



in the wall, we obtain 



39 



for the solid wall. Note the introduction of h in (19); see remarks 



following (17). 

b. Fluid; Initial and Boundary Conditions 

(1) Initial Condition t = 0. 

0* = 0 0 ^ x < oo , t = 0 

(2) Boundary Conditions at x = 0. 

0=1 x = 0 , t > 0 

(3) Boundary Conditions as x — >co . 

0 > 0 x — v oo , t > 0 

c. Solid Wall; Initial and Boundary Conditions. 

(1) Initial Condition, t = 0. 

T = 0 0 ^ y 2 ^ b ’ t = 0 



(2) Boundary Condition at ^2 = b. 




(3) Boundary Condition at ^2 = 0. 

- K = h (0"- T) y = 0 , t> 0 

m 2 

dy 2 



This system which is similar to that obtained in subsection IIDld is 
somewhat easier to deal with. Accordingly, it is treated first, in VA 
while that for the system IIDld is dealt with in VB. 



III. INTEGRAL EQUATION FORMULATION 



A. INTRODUCTION AND PLAN OF ANALYSIS 

Equations (II - 3 ) and (II-4), with K£ a = Kfj. = K ^ and K ma = K m ^ = K m 
so as to provide isotropic thermal conductivity, involve two dependent 
variables, 0(x, y, t) and T(x, y, t), each involving three independent 
variables, along with a formidable number of boundary and initial 
conditions. Similar remarks hold for the case of cylindrical geometry. 
These are complicated systems to deal with and we have achieved 
solutions which are useful for calculating purposes only by making 
various simplifying assumptions. These solutions are given in other 
sections of this thesis. In this section we will represent the system 
of equations (II-3) and (II-4) in another form which involves only one 
unknown function of two independent variables. No additional conditions 
are required for the complete specification of this unknown function. 

The form of this equation, (70), is that of an integral equation of what 
is known as "first type" (that is, the unknown function does not appear 
outside an integral sign) with complicated kernels and forcing functions. 

Although the method of procedure is essentially the same for the 
case of cylindrical geometry as for the case of slab geometry, the 
details are sufficiently different to justify treating the two cases 
separately. We treat the case of slab geometry first and in subsection 
HID we return to the case of cylindrical geometry. In each case the 
method of procedure is essentially as described in the next paragraph. 
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We define an unknown function (x, t) which represents the 
(unknown) temperature of the fluid at the interface and a similar 
function y(x, t) which represents the unknown temperature of the 
solid at the interface. For the fluid we solve the problem of finding 
the interior temperature 0(x,y,t) when the surface temperature is 
given by ^(x, t), subject to other appropriate boundary and initial 
conditions. It is convenient to consider this solution as the sum of 
two solutions, one involving a homogeneous boundary condition at the 
interface and accommodating the boundary condition at the end, x = 0, 
and a second solution involving a homogeneous boundary condition 
for x = 0 and the unknown function <j>(x, t) at the interface. The 
solution to the latter problem is expressed in terms of an integral 
involving the function^(x, t); it is obtained by use of a procedure 
which is due to A. N. Lowan [20, 21]. We proceed in the same way 
to express the solution for the solid temperature T(x, y,t) in terms 
of the unknown surface temperature ^(x, t). The interface heat 
balance condition provides two relations, one of which may be used 
to express the function *l|/ (x, t) in terms of the function <f> (x,t). 

The other relation then becomes the integral equation for the unknown 
function t (x, t). 

B. AUXILIARY PROBLEMS (SLAB GEOMETRY) 

We will deal with the problem of the temperature distribution 
in the fluid first. Then, by taking the flow velocity as zero and by 
using subscript m (for the solid) rather than f (for the fluid), and 
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by changing the designation of the surface temperature distribution 
from $ (x, t) to f(x, t), we can very quickly obtain the corresponding 
solution for the temperature distribution in the solid. We immediately 
write 

0(x, y, t) = u^x, y, t) + u 2 (x, y, t) (1) 

where Uj(x,y,t) are the solutions to the two auxiliary problems 
considered in the remainder of this subsection. 

1. Fluid; Homogeneous Boundary Condition at Interface 

For simplicity in this subsection we will omit subscripts, 
and will replace them when we combine solutions: cX will denote 
0^£, u will denote uj, y will denote y^ (Cf. Figure lb). Thus the 
problem of interest is 

Hr . ~£>ia> V Da / 2 \ 

bx 2 - -£ X bx ~ « bt 

The boundary and initial conditions are 



u = 1 


x = 0 


0 Cy^a 


t> 0 


(3) 


u — 0 


x -> oo 


Oly 4 a 


t|- 0 


(4) 


o 

II 

2 


O^x < oo 


II 

o 


t> 0 


(5) 


"^ = 0 


0*X < OO 


y = a 


t? 0 


(6) 


u = 0 


8 

V 

X 

II N/ll 

o 


o<.y= a 


t = 0 


(7) 
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Using an overbar to denote Laplace transform with respect to time, 



we obtain 



J- - 

^ u Su. 
+ bx* 






_V. 

* bx 



-i-a 



( 8 ) 



This incorporates the initial condition. All the boundary conditions 
will be the same in Laplace transform form except that at x = 0 

which becomes \Ji = Thus assuming a solution of the form 

S 



— -A* 

U = A e " sin(|jy) 



(9) 



we observe that all conditions are satisfied except that for x = 0, 
and the condition of zero slope at y = a. Imposing the latter yields 
the characteristic values 

jU n = (n+1/2) IT /a n = 0,1,2,... (10) 

Because of linearity we can write 



co . 

An e" ***(/U) 



( 11 ) 



n*o 



Substitution of (11) into (8) yields 



2<* 



( 12 ) 



The negative root is discarded in order to assure boundedness as 
x — v oo • The condition at x = 0 and use of the orthogonality property 
of trignometric functions yields 
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A = — 

« s 



/V 



( 13 ) 



Substituting from (13) and (12) into (11) we obtain 



oo , 

<*p [js - * yf* + 1- 

«=o 



(14) 



let 2 



V 



where 

A7. M r^' (15) 

The inversion of (14) may be obtained using standard tables 
[29]. Thus 

Co 

U(x.y,t)= Y^y) X^x, t) (16) 



where 










(17) 



V* 



~/ i - 

€ 



• e r/c (' 



x-z«/V 

l^Zt 



* 

/i* 






(18) 



In these expressions, although there is explicit dependence upon the 
integer n, this is not reflected in the notation on the left. This is 
done merely for simplicity of notation and will also be done with 
similar definitions to follow. In every case the symbols are used 
under a sign of summation with respect to n. 
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For future purposes we evaluate the partial derivative of 
(16) with respect to y and evaluate at y = 0. 

oo 




where 

y'(0) = 1 (20) 

As the reader may come to appreciate, when the many preliminary 
results are assembled into the final equation, a compact notation is 
convenient. Functions and and constants jJ ^ are the first of 
several such notations which will be introduced. 

2. Fluid; Non-homogeneous Boundary Condition at Interface 
In what follows U stands forU^(x,y,t) and it is governed by 
the same equation (2). However, the following changes are made in 
the boundary conditions (3) and (5) 

U = 0 (rather than u = 1) x = 0 0<y<a t > 0 

u= t(x, t) (rather than u = 0) 0 ^.x^oo y=0 t^O 

where | <p (x, t) | =M<<& for all x and t. We further make the 

following substitution in order to drop out the convection term. 

u(x, y, t) = W(x, y, t) exp « Vj: ) (21) 

U°<. tfoc J 

Then equation (2) and the appropriate boundary and initial conditions 
reduce to the following equation and conditions for W(x,y, t). 
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( 22 ) 



cfvy , tfw _ j_ W 
bx 2 « bt 



w = o 

W<oo 

W = 0(K,tj 

** 

W- o 



x=o 0 < 2 <a J 

X— >oo j o < a. j 

o^y < oo y ^ = <x ? 
0 <X <oo J 0 <y < 0 . J 



i>o 

t l° 

b t° 

t>0 

i : =o 



where 



3(X,tJ = <h<p 



- VX + vV 

2* ^oc 



(23) 

(24) 

(25) 

(26) 

(27) 

(28) 



Again taking Laplace transformations with respect to t and accounting 
for the initial condition we obtain 



-vt— • O.- _ 

dw . ow - 1 vV 
*x 2 * 



(29) 



All the boundary conditions will be exactly the same as above except 
(25), which becomes W = g (x, s). 
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Following the method of solution used by Lowan [20, 21] we 



consider the function 



OO 00 



O 0 



(30) 



where 



F (a,u;U [X. ( a-8)J 

* ' cost (a a; 



(31) 



and 

X 2 = yu 2 + ~ (32) 

Note that j* and ^ are dummy variables of integration. Direct 
substitution shows that (29) is satisfied. It is easy to verify that all 
the boundary and initial conditions are satisfied. In particular we 
note that 

CO CD 
0 0 

by the Fourier identity. As regards boundedness 

OO OO 

< 2p. Ava px. aL*. p'g djj = liU <oo 

This completes the proof that (30) is the desired solution of (29) 
which satisfies all the auxiliary conditions. 




48 



Continuing to follow the procedure described by Lowan [20,21], 



we define 



0 " 



(33) 



i. e. , (a,y;t) is the inverse of F (a,y; A ). The subscript jj. 

indicates the dependence upon the parameter y- ; the dependence of F 
upon s is made clear by expressing A in terms of s, as above. The 
convolution theorem then gives 

oooo t 

°J|A (34) 

Q O q 



Further defining 




(35) 



then 



7>t 



- Ff, (M'*) 



Thus substituting from (36) in (34) yields 



(36) 
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OO CO t 

f ^l a J ft (^ t- >1 jc^ 

o c> "L ~^ r ^ 



(37) 



If we now take the Laplace transform of (35), we obtain 

. i YU) 



s £(s; 



where 



(38) 



Y(sj = [x(a-^)J 



Z is) = CtfS-& ( Xa) 



(39) 



The singularities of (38) are simple poles corresponding to 

the zeros of cosh( X a )*> of course, these are functions of the variable 

jX . Note that cosh (A a) = 1 + V^/2 ! + X^a^/4! + ... which converges 

everywhere in the finite plane. Thus, although X itself involves a 

2 

square root , X does not, and there is no difficulty to be en- 
countered dealing with branches. Thus, the residue theorem affords 
the desired inverse of (3 8) 



zS + 1 



QO Sn'f’ 

Y(q) A V 1 -€. 



Ylzo 



m. 



(40) 
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The summation in (40) extends over the roots of 

cosh ( A a) = 0 (41 ) 

namely the numbers 



- P +1 * (n+ ,/z ) 7T = 



n =. 



1,2, 



(42) 



(Do not confuse the constants f*n of (10) with the continuous variable 
f of (30)). Thus 

S„ B _o< [(*+. i/ z ) + fj 1 -! 

L r J (43) 

= - «< (/'a + /“ J 

where S n are the poles of (38) in excess of the simple pole at S = 0. 
Note that they are functions of the variable ^ . Therefore the formal 
inverse of (38) is 

Cos ^(HaJ ^ L ( + p % ) 

n=o 



$ 






(44) 



The following derivative is needed in (37) 



"b-t 



Co 

it c< * i^'^) 



n^o 



-*(ji ■*!?)* 

e 



(45) 
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Hence, substituting from (45) in (37) we obtain 




In order to complete the analysis, substitute (28) in (46) and then 
substitute the result into (21). Thus we obtain 



( 46 ) 



*'*>*>*> 




o o 



f(^H) 



/ 



co 2 

i , 

e /iZ* Aft 



) 

j<M§ 



(47) 



We are able to carry out the integration with respect to JJ so that we 
finally obtain 






» oo t 

2^Y t (sJ jy ■Kf 



ri=o 



o o 



(48) 



where 



YJa)= ^(-o n cos 

u ( , , j"5T pi &< “t (%U-‘- o( t J ») 

K| e . e 



1 


r 






'-Os + xfl 


ft ' 


e^.p 


r J 


— e-txp 


J > 

4 



(49) 



( 50 ) 
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which is evaluated below. 



For later purposes we need 

OO CO -t 






3=‘ 



-dLL 



| = Y , (*;? 

a=° fe, 0 ° 



(5i: 



In order to assure that there is no misinterpretation, we 
remark that the integration with respect to ^ is from zero to infinity 
and that with respect to extends from zero to t. 

3. Combination of Solutions for the Fluid 
Simple addition of (16) and (48) yields 



* 00 * 

0(x, 3 ,*) =V t ff' 

nso 



(52) 



Likewise by superimposing (19) and (51) we obtain, 



( 



f^] =^[x,(*>‘)+//'K?,t-i)K,(*,! (53) 

3. 0 ° ^ 



*n=o 



o o 



where, having restored the appropriate subscripts. 



Y,<a,) = c ^ h - 






Vx 
I Zo< f 

— e f 



* 

e". eric ( ) + e • erfc ( ^ tUd) 

' 2 Y V / ' 



r" 






(54) 



(55) 



( 56 ) 
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( 57 ) 



Ki ( * j ^ t j 




i 





\ 






f 


_ l 




, i “ 


\ 






- exp 






< e x p 




^°yl 
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Also, of course, (10) and (15) define jj and . 

It is well to pause at this point to review the physical nature 
of the problem for which (52) provides the solution. This is given 
below. The differential equation is 



V dj? _ j_ 

lx 7 -' 7>x + dt 



The auxiliary conditions are 



CO 

X 


X = O 




0 i i a 


S 


t > 0 


0 < CO 


X-*- oo 






J 


t > 0 


e * 


0 : 4 : X < 00 




c/,= ° 


> 


t> o 


a 3 


oi X < 0 o 




3,-* 


3 


t> o 




0 < X < oo 


) 




2 


o 


Solution for the Solid Wall 











The solution for the similar problem for the solid wall is 
easily obtained from the preceding solution for the fluid. All that is 
necessary to do is to write T in place of 9, Y i n place of , replace 
subscript f by subscript m, replace a by b, take V = 0, and 
introduce a new symbol (defined below) in place of . 
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oot 

T(*jz> t ) = ) y 3 (*) Xz(^t) + f [rdf m)x 2 jj-i j * 

H-oi 



(58) 



where 



= 



/^W 



(59) 



X 2 ( x -> *■ ) = -r 



^ ^ c ( 7 $rke'' e V ( W^i 



(60) 



Y«W = 



_ (-0 CoS [x«(t-a)] 






(61) 
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(62) 



An = (n + 1 / 2 . ) -£• 



(Do not confuse X n above with X A used temporarily in (42)) 



(63) 
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Op CO t 

3 r° ^ — ' L 0 c 

> i=o 



(64) 
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C. DERIVATION OF INTEGRAL EQUATION (SLAB GEOMETRY) 
Equations (52) and (58) thereby give the fluid and pipe wall 
temperature distributions respectively. However these solutions are 
in terms of the arbitrary interface temperature functions 4* (x, t) 
and Y (x, t). We now return to the actual boundary conditions that 
prevail at the interface, viz. 

= -R (^-y) jj = ° > <o ° > t ^° (65) 



and 



. i (t-y) 



fix.- 0 j ( 66 ) 



(Note: The signs on the left hand side of (65) and (66) are consistent 

with the choice of coordinates used in Figure lb). Solving for Y from 
(65) in terms of , we obtain 






(67) 



Substitution from (53) in (67) yields 






it 



B 



°? t 



0 O 







( 68 ) 
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Furthermore, substitution of (68) in (66) yields 




co t 

o o 




( 69 ) 



Substitution of y' (x,t) from (68) into (64) and the result into (69) 
yields 



oo 00 r oo 

n=o 0 ° m~& 



o o 



^ j <{>(«; *-1 -A) 



o o 



+ 



(70) 



CO r CO V 

0 0 



= 0 



It should be clear that we could, with equal facility, express the 
problem in terms of y (x, t) instead of (^> (x, t) as above. The 
distinction is trivial. 

Equation (70) is the desired integral equation, the only unknown 
in which is the function <f> (x, t) which represents the surface tem- 
perature of the fluid. 

Equation (70) is considerably simplified if there is no surface 
resistance. Letting h, the heat transfer coefficient, approach 
infinity, reduces (70) to 
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Oo I oo t 

o-£ |X|(x,t) + + J 

n=o [ f 

[Xi( x >?-» ) + "£~ x 2 (>cf , ^1) 



A d )j , 



( 71 ) 



for h = co . 

D. AUXILIARY PROBLEMS (CYLINDRICAL GEOMETRY) 

For the case of cylindrical geometry, the treatment is quite 
similar to that which has preceded. However, it is not possible to 
derive the equations for the pipe wall by simple manipulation of the 
equation for the fluid as was done above. Thus it is appropriate to 
consider four auxiliary problems similar to those dealt with before, 
two for the fluid, arriving at equation (165), and two for the solid, 
arriving at equation (170). These equations will then finally be 
combined so as to arrive at the desired integral equation. By suitably 
redefining the several functions which are involved, it will be possible 
to express this equation, (181) exactly in the form of (70). 

1 • Fluid; Homogeneous Boundary Condition at Interface 
Consider the problem mathematically stated as 

v_ _ i 'du. 

bx z * dx r 2r dr z * bt 1 ' 
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with the boundary and initial conditions: 










U = 1 X= o 


j a<r <a 




-t > o 


(73) 


U. s o X — > <*> 


j »< c <a 




t>D 


(74) 


^ = 0 a <x < ao 

br 


) V -o 




t><5 

r 


(75) 


U.=* O o<x < oo 


•> <•?. a 




t>o 


(76) 


U - o o<X < co 


o<r < <k 




tx TO 


(77) 


Assuming a solution of the form 










u = A X,(x 


>*) 






(78) 


and substituting (78) in (72) yields 

tf-Xi _ ax. _ = 

Dx* * 


_L dXi 

* bt 






(79) 



Equation (78) already satisfies condition (75). Condition (76) yields 
the characteristic equation 

«!( P Gj = 0 (80) 



the roots of which are denoted by n=l, 2, . . . . Because of the 

linearity we can write 

00 

u= ^A n J 0 X) (81) 

The solution of (79) is obtained by taking Laplace transform with 
respect to time, which yields a second order homogeneous ordinary 
differential equation the solution of which is 
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(82) 



X,(x,t) = B -**p 






In (82) the exponential with positive exponent has been discarded so 
that condition (74) may be satisfied. Substitution of (82) in (81) and 
absorbing B into A n , we obtain 



oo 



g-xfef+A 1 * 5'' 



(83) 



Thus at this point we observe that the inverse of (83) satisfies (72), 

(74), (75), (76) and (77). Condition (73) provides a means of evaluating 
A n . Taking x = 0 (73) gives 



oo 



. 1 - = ^ A n Jo ( 

n- 1 



(84) 



multiplying by r J 0 ( /? P )dr and using the orthogonality relation of 



Bessel functions; 

a 



fjAP>nr)J 0 (P n r)rJr = ± J, 



(85) 



if n=m and zero if m, one easily finds 



f\ _ I 2 !_ 



( 86 ) 



Thus we now have 




( 87 ) 
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The inversion of (87) is tabulated [ 29 ]; viz. 



a (r,x,t) = 



where 



I R,(r) 



X, (x,t) 



( 88 ) 



R,(r)= _L. '1 
P" JAt-o-) 



= i r 

<x 



e • erfc ( X T0) + e'* er t c ( % '0^) 



(89) 



90) 



<*£ »/?. + ( v / 2 « / 



For future purposes we need the partial derivative of (8 8) 
with respect to radius evaluated at r = a. 



-bUL 

3r 



r= a 




*=i 



(91) 



2. Fluid; Non-homogeneous B.C. at Interface 

The problem considered here is the same as that in the 
preceding subsection except that the following changes are made to 
the boundary conditions (73) and (76). 

U = 0 (rather than u = 1) x=o J 0 < r < a f t>o 

a = 4 ) (x,t) (rather than u = 0) o< x<Q ,, r,A j t> o 

where | <p (x, t) | <M < 00 for all x, t. A substitution of the form 

Xjt) ~ 7^ ) (92) 
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reduces (72) to 



cfw + x + t?v - i w 
^x 2 - r dr ^r 1 - «• &* 



and all the boundary and initial conditions will remain the same 



except that at r = a we have 



W(a,x, t) = <$> (x, t) exp + ) = g(x, t) (94) 

Taking Laplace transform of (93), we obtain 



5 w + X ^ X = X w 

dx* r <jr dr 2 - * 



(95) 



Again all the boundary conditions will be the same except (94) which 
in Laplace transforms will become 



W (a ^ (x*s) 



o <oo r= a (96) 



Following the method of solution used by Lowan we consider 
the function 



oo 

^ffF( a ' r ' X ) 3 (97) 

o o 



where 



and 



F (<Xr;X) = LiAf. } 

X>(xaj 




(98) 



(99) 
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Note that y> and ^ are dummy variables of integration. Direct 
substitution shows that (95) is satisfied as well as all the boundary 
conditions. In particular once again we note that 

co oo 



^ S )- ^ ff £ ^ 5 j A 

0 o 



( 100 ) 



= 3 (*;!) 



by the Fourier identity. As regards boundedness, 

oo oo 

|w(r, X;S) | ^ ^ J j j. X X Aav\ jjL'fZ d ju d^ = m C oo 

o o 



( 101 ) 



This completes the proof that (97) is a solution of (95) which satisfies 
all auxiliary conditions. Let us now define 
co 



Je 1 F r K'V t-JJf = F (*, <•)*) = F fa*-; f^f) 



( 102 ) 



i. e. , (a, r;t) is the inverse of F (a, r; X ). Note again that the 

subscript ^ indicates the dependence upon the parameter ^ ; the 

dependence of F upon s is clearly indicated by expressing A in 
terms of s, as shown in (102). The convolution theorem then gives 



OO co l 

W(r,x,*J = f (103) 



o o 
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Further, defining 



t 

Gfi (ajrjt) = Jffj (a. J rjir)Jt 



gives 



= Fu (a^rjt ) 
at y 



Substituting from (105) in (103) yields 



O 0 0 



If we now take the Laplace transform of (104), we obtain 



oo 

£. < y(<*>r';t)|. J r; tj Ji 



- J- f(“^<*;X) 



= i Y(Q 
s 2(sJ 

where 



Y(*)= M*r) ; 2 (s)= I«(A<x) 



( 104 ) 



(105) 



(106) 



(107) 



(108) 
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Using the (temporary) substitution \ - L fi 9 (98) may be written 



f - Mill 

J. 



( 109 ) 



However J 0 (^ ) is an entire function having only real zeros. Thus 
the only singularities of F are poles at the real roots, f? n of 

J oip*) = 0 (110) 

The corresponding values of s are: 

s (A pS) < m > 

Following Lowan, we write the desired inverse as 






Y(°J 

2(o) 



+ 




Y(s-; 



s*t 

e. 




= s 






. oo » 

Io(D _ 2 y p*J,Me 

[far') j, {p**) 

»=i 



( 112 ) 



Thus taking the derivative of (112) and substituting into ( 1 0 6 ) yields. 



W<r,x,t)« -i£L 







(113) 



It is possible to carry out the integration with respect to p. 
in (113) and obtain 
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( 114 ) 



00 rr 

U(r,x,t) = VV(r)JJ 



where 



R 2 (r) = 



JJM) 






z -t(& + «/£) 



* : t : _L . i* 

W * ft 



■ e 



1 e-,C P 



4* »\ 



Hfe+jO 






(115) 



(116) 



Note that K^(x,^ , t, *\ ) given in (116) differs from K|(x,^ , t, \ ) given 
by (50) only in the eigenvalues; otherwise they are identical. 

For later purposes we need the partial derivative of (114) 
with respect to r t evaluated at r *» a* 



bu 

dr 



OO oo t 



(117) 

«=/ ° * 



3. Solid Pipe; Homogeneous Boundary Condition at Interface 
Consider the problem stated mathematically as 



d^q f i ^ ( r ^_u \ _ ± 7)0. 
dx 2 r hr'' 7>r * <* <)t 



( 118 ) 



with the boundary and initial conditions: 



66 



u= i 



U- - o 






u= o 



a = o 



x=o , 4 <r ) t > 0 ( 119 ) 

x~+co j a<r<t , t>o (120) 

o^x <00 J r= b > t> o (121) 

0 < X<00 J r= O- > t-Zo (122) 

o <oo j> atr <6 j "t-o (123) 



We first define the following notational device for a linear combination 
of Bessel functions of the first and second kind of an integer order. 

Let 

b-i) = J K (^)y n (vt>)-Y K (yr)J rt (!ri)j (124) 

Nowassume a solution of the form 



u - A (r, bj x) 'XiixA) 



Substitution of (125) in (118) yields 




« J-t 



The function (125) already satisfies ( 12 1), i. e. , 



(125) 



(126) 
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= o 



(127) 



t>r I 

T=b 



* £ ( b, b; * ) 



Condition (122) yields the characteristic equation 
a&i (*, b* ) =0 (128) 



the roots of which are denoted by Y n , n = 1, 2, . . . . Because of 
linearity we can therefore write 



oo 



a = 






(129) 



To solve (126) we take Laplace transform with respect to 
time. This yields a second order ordinary homogeneous differential 
equation, the solution of which is 



X 2 (*; s ) = 



8 -^xp 




(130) 



In (130) the exponential with positive exponent has been discarded so 
that condition (120) may be satisfied. Substitution of (130) in (129) 

and absorbing B into A^, we obtain 

oo 

e € l (x';bj'x) e*p 

Thus at this point we observe that the inverse of (131) satisfies (118), 
(120), (121), (122) and (123). Condition (119) provides a means of 
evaluating Taking x = 0 (119) gives 






(131) 
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(132) 



OQ 

s = ^ ^ ^ ^ 



n -I 



It is not too difficult to show that (124) has the usual orthogonality 
properties of Bessel functions with the weighting function r. Thus 
we merely present the result as 

/ » y " ' ' - (133) 

i a & W.)- 

if n = m and zero if n ^ m. Use of this orthogonality property yields 

& (a. hi r.) 



A = J. 
" s 






z 



(134) 



Substitution of (134) in (131) yields 



U 



-Z 

Tl=r| 






j ^p(-x y^+r ) 



(135) 



Note that the parameter s appears only in the last term. The inverse 



of (135) is tabulated. Hence, 

oo 



u 



(r,x, t) =^~R 3 (r) X^(x, t) 



(136) 



where 






R ( r ) = J2L [fl* kj SVi ) o^t b; ]fn j 



(137) 



X.Cx.-tj = e^Vfe + e'\ rf c (138) 
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The partial derivative of (136) with respect to r and evaluated 



at r = a is 



du 
dr . 



r=£ 




X z t) 



where 



(139) 



M*) = 



a 



v* 



(140) 



This result will also be used later. 

4. Solid Wall; Non-homogeneous Boundary Condition at Interface 
We consider the immediately preceding problem with the 
following changes in boundary conditions 



U = 0 


(rather than u = 1) x = 0 


a< r* < 6 




u = Y (X, 


t) (rather than u = 0) x - 0 


r =• o- 


J 


where | \j/ ( X , t) J 


(oo for all X and t. After Laplace 


transforma - 


tion we have 









, ± i ( r 15 \ 
r 2 )r ' dr / 




(141) 



with the conditions 



U = o 


X =o 


a < b 


(142) 


U ( oo 


y — v oo 


a <r < t 


(143) 
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= o ° < X <03 , r=b 

br ' 

U = Y(x,sj a ^ x < & j r ~ a 



(144) 

(145) 



The solution of (141) may be constructed following the procedure used 

by Lowan. Thus, consider the function 

oo oo 

“ - \ ffei (i46) 

0 O 



where 



F( n Av;\) = 



r„(*r) Ki(>i>) t K 0 (Ar) r,(M>) 

I.(H k,(Ab) + K„ (\a) I,(AbJ 



(147) 



and 



^ p + s/« 



(148) 



(I and K are modified Bessel functions of first and second kind of 
order n). Observe that jJ and ^ are dummy variables of integration. 

Direct substitution shows that (141) is satisfied. It is also 
easy to verify that all the boundary conditions are satisfied. In 
particular we note that 



oo oo 

U = A J'p y ^ ^ 

O O 



s Y(x^s) 



( 149 ) 
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by the Fourier identity. As regards boundedness, 



co oo 

A ^'^ x if ^ " (150) 

o o 



This completes the proof that (146) is the desired solution of (141) and 
satisfies all the auxiliary conditions. Again we follow the procedure 
indicated by Lowan, and the sequence of analysis is precisely that 
employed in subsection IIIB2 for the case of slab geometry; the present 
case is different only in details. Accordingly, we abbreviate the 
treatment and write immediately, 

oo co t 

0 o 0 L 



and 



= / 1 fi F ( a, b , = /' 



l Yi£ll 

s 2(0 J 



(152) 



where 

Y(S) = J o (Xr) K,(xb) + K*(Xr) r,( Xk) (153) 

2(0 * I 0 (xajK,(Ak) + K 0 (Xa) I, (M>J (154) 

Using the temporary substitution X= iY , (147) may be written 

F( a a,*). iMW-UnKM _ Ajf^sL ( i55) 
J,MY,(n)-j,(n ) xM A 
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where we have used the definition (124). It is reasonable to assume 
that the function (a,b; H ) considered as a function of t has 

only real zeros. Thus the only singularities of F are poles at the 
real roots, y , of 



6 j TfnJ 



(156) 



The corresponding values of s are. 

S n =- (** + f* 1 ) (157) 

Following the procedure used by Lowan, for the inversion of (152), 
we obtain 



= "M + Y- 

H(»; L, 



CO 

Y(s« ) € 



sf5Jl 

1 ^ 



(158) 






G r (a, b, r; tJ . Io(HdK«(Ht>)-fKo(pr)X,(pO + 

L(P a )K,(p*>) + K.(pa) I, (pa) 



CO 



m„) 



V 3 l_ 

tf+p 2 - 



(159) 



M _ t (^jbjYn) + a ,b( (0>i>.>K*)j 

Thus taking the derivative of (159) and substituting into (151) yields. 



co L 

K(tx ^* Twfl J ) K2.(Xj^ , tj<^) (160) 



where 
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(161) 



R 4 (r) = 



Xn (fj bj Yn ) 



-c<r M i 



Kz(*^mj= 



exp 



-Jtd 



- «-*p 



'-M 1 



(162) 



For future use we require the following 



w 



oo , 00 x 

= V R 4 ( *) /JV C ^ * - n. ) K 2. (* ^ , t, a) d ^ 5 

^ 0 0 



I 



(163) 



where 



RqH = 



0^0 6j y»i) 



(164) 



-t a 



5 . Combination of Solutions (Cylindrical Geometry) 
a. For the Fluid 



Simple superposition of the solutions (88) and (114) yields 

" cot 

e(r,x,t>2^ ta(dX,(x,t} + Rl<.r)ffo(hH)K^,t,' l )d'i<)?- (165) 

n=i 00 J 

where having restored the appropriate subscripts 
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R,(r) =R t (r) = j 

pn J, (M) 



( 166 ) 



x,(x,tj= ± h 

& 






i . *> 



,, , K7 « ! £S>-’5) <*V~A«() 

K t (x^,M)-#j==- e f • e 



• < 



exp 






Z 



— ^Kp 



-0[W 



2 1 



Similarly superposing (91) and (117) we obtain 

<*> a> t 



' be 



= -^T [x,(x,t; +J JtO^t-a) 
r=a *7! 00 



b. For the Solid Pipe 

Again simple superposition of (136) and (160) yields 
op a? “t 



V|=r« 



O 0 



where again restoring the appropriate subscripts 

Rj( r ) - A ^ bj *” ^ 






Rq(fJ=. o^l ( r -» kj W 






i -* m *„*l 

1 / / I V 'Q 

Kz(M^^ = >nr 









-(T-xJ 

« v ^ 



- ^ p 



-M 

L ^"A 



xy 



(167) 

(168) 



(169) 



(170) 



(171) 



(172) 



(173) 
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similarly superposing (139) and (163) yields 



CO co t 

+ (174) 

o r J r* a ^ ^ o ^ > 

n=l 



where 



/ 



R,(aj 



R,(«) = 



1 



-1 



Jb #^i (bj t>3 ifn ) 
a «£,(»,bjQ 

Yn Z 

b 0 ^>( g -’bjy n ) ^ /? 



(175) 



(176) 



E. DERIVATION OF INTEGRAL EQUATION (CYLINDRICAL GEOMETRY) 
We now proceed as we did in section IIIC for the case of slab 
geometry. Equations (165) and (170) give the fluid and solid temperature 
distribution each in terms of the arbitrarily preassigned function 
cj> (x, t) and y (x, t). If we now return to the actual boundary condition, 
namely that of continuity of heat flux at the interface, we have 




2© 

ir* 






-K 



/ r=a 






(177) 



and 




21 



= -R (?-y) 

'Vo a 



(178) 
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and solving (177) for ^ in terms of cf> , we obtain. 



y =. 

£ 



£3 

'br 



+ <j> 



( 179 } 



r„a 



Substitution of (169) into (179) yields. 

00 CO t 



Y= f - jj (180) 

n* 1 o o 

Further substitution for^ from (180) above into (174) and the result 
into (178) yields 

oo f co t 

KyyA < $1 ( &) O + Rq (&) jtj r^'j • 

^ O n 



= i v 






o o 

OO t-*t 



(181) 



OQ 00 t L 

> r»=i 0 0 

oo 



’Kj (t >Pj *" y i> *) /+(W-i)Ki(*^M)V5]=<’ 

JJ rjssl o * 



The same remarks as those following equation (70) for the integral 
equation of slab geometry may be made here too. Similarly therefore 
if we let h approach infinity, we obtain 





Ra(«) Kl(X;f; M) 



/ = ^ 



(182) 
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This is equivalent to (71). Note that (182) can be written more 
compactly as 



n=( o o 



o 



( 183 ) 



where 

X ( x jt) = —■ - $3 ( 0-) X2.( x -»^) + X((x^'t > ) (184) 

K i 

K K,(x,^M) + ^ (185) 

A similar abbreviation may be made for (71). 

The reader should note that (70) and (181) are of exactly the 
same form; their only difference is in the particular form assumed by 
the kernel and driving function and in the fact that the counterparts of 
Rg(a) and R^(a) appearing here are simply unity in the case of (70). 

F. REMARKS ON SOLVING THE INTEGRAL EQUATION 

It is obvious that the integral equation (70) or (71) for the case of 
slab geometry, or (181) or (182) for the case of cylindrical geometry 
will not be solved easily. It involves multiple integrals, one with 
variable upper limit and one with fixed infinite upper limit. The kernel 
and driving functions, although well behaved, are quite complicated. 
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One of the ways of dealing with such a complicated problem is to do 
what has been done in this thesis, namely, make additional assumptions 
so as to simplify the problem to the point that it becomes analytically 
tractable. Our attention has been exclusively directed in this way and 
we have made no serious attempt to deal with these fundamental integral 
equations. However, the following remarks and conjectures may be of 
interest. 

All the analysis herein shows that while <f> (x, t) is a smooth 
function, almost all its variation takes place in the immediate 
neighborhood of the wave front x = Vt. It is conjectured that a contour 
map of (x,t) would be similar to that shown in Figure 2. In this 
figure the contours are shown as straight lines merely because the 




Figure 2. Conjectured Contours of <f> (x, t). 
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nature and degree of curvature are unknown. It is tempting to consider 
replacing the boundary and initial conditions = 1 along x = 0 and 
^ = 0 along t = 0, by the conditions $ = 1 along OA and = 0 along 
OB so as to provide a field in which the "action" of <J> (x, t) is more 
evenly distributed. This might perhaps be useful in an attempt at 
solution via Gallerkin’s (or a similar) method in which a solution is 
assumed which involves unknown scalar constants and these are 
determined so as to minimize an error functional. In such a procedure, 
it is evident that simplifications would result from assuming the con- 
tour lines to be straight. 
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IV. METHOD OF PERTURBATION 



In this section we employ the method of perturbation to treat the 
problem described by equations (II- 10) and (II- 11) and the associated 
boundary and initial conditions. The reader will recall that these 
equations resulted from having taken radial (or transverse) averages 
of temperatures in the fluid and in the pipe wall. Thus radial (or 
transverse) conductivity is, in effect, absorbed into the surface heat 
transfer coefficient h* defined in equation (II-9). The influence of 
axial conductivity appears explicitly in the following analysis. 

An outline of what will be accomplished in this section may be 
helpful to the reader. First, an ordinary perturbation expansion is 
made in terms of the parameter 6 = 1/V. It turns out that only even 
powers of £ are involved. The zeroth order solutions are identical 
to well known solutions described in section 15-3 of [7]. The effect 
of longitudinal conductivity does not appear in this zero order solution 
but it does appear in subsequent solutions. It is clear from the sequence 
of higher order perturbations that the solutions fail at and ahead of 
the wave front, x =Vt. Severe singularities are encountered at the 
wave front and the solutions are identically zero ahead of the wave 
front. 

Thus a singular perturbation development is introduced to explain 
the behavior ahead of the front. The lowest order solution for this 
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singular perturbation expansion is similar to the known solution for a 
moving fluid, with diffusion and with an isothermal boundary condition. 
Terms in this singular perturbation development involve unknown 
functions of time which are determined by a process of matching with 
the solutions valid behind the front. 

It should be remarked that the system of equations dealt with in 
this section is obtainable from either slab geometry or cylindrical 
geometry. The process of radial or transverse averaging of tem- 
peratures erases the distinction between the two geometries. 

A. ORDINARY PERTURBATION ANALYSIS 

Equation (II- 10) and (II- 11) together with associated conditions 
are repeated here for reference purposes. The notation has been 
slightly changed; we omit the asterisks which indicated radial averages 
and we introduce the convenient constants b^ and b^ defined below. 

In the fluid, we have 




( 1 ) 



where 




= if 



( 2 ) 




subject to the conditions 
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9=1 


x = o , 


H- 

UV 

o 


0 — o 


X — *■ co j 


t ^ O 


6 -o 


o< X<oo j 


£ - o 



For the solid we have 

« £l « ~£ _ bi(e-T) 

rn bt 

where 




subject to the conditions. 



( 4 ) 



( 5 ) 



T= i 
r-o 
T= o 



x = a ^ t >o 

X — ► oo j 't > o 



o <X y t-o 



( 6 ) 



It is convenient to normalize with the following change of 
variables 
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t a t 



( 7 ) 








1 = x /y = 


£ X 


; t 


= l/y 




o that equation 


/ 

(1) and conditions 


(3) become 








2 t \q 

o< € — ~ 

f 


tf 


= 5f + b,(e- 

~bt \ 


■T) 


(8) 




9-1 




V s0 j 


t^O 




< 


9—^o 




*-00 .> 


t >0 


(9) 




9 = o 

j 

\ 


( 


? < oo 


t = 0 




.nd similarly equation (4) and conditions (6) become 






u itl 


_ ~bT 


- b 2 (e-T) 




(10) 




T-i 




1 


tz o 




< 


T-v o 




*{-+00 


t> O 


(11) 




T=o 




0 < oo 


t — 0 





Equations (8) and (10) with conditions (9) and (11) can be approached by 
the use of Laplace Transformations as was done by Hiep [13], but 
there seems to be no way of solving the transformed equations or 
performing the inversion other than by strictly numerical procedures. 

/ 
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This system of equations may be considered as depending on the 

2 

parameter £ which appears in (8) and (10). For large values of V, 

this parameter is small so that it is reasonable to consider a per- 

2 

turbation expansion in £ . 

2 

(Note that £ itself does not appear in (8) through (11), only £ ; 

2 

thus the expansion we develop is in powers of £ , that is, only even 
powers of £ are involved. There are other ways of showing that odd 
powers should not appear in this perturbation expansion, but the one 
given here is sufficient. ) Thus, 



e) =T 0 (>i,t) + + (13) 

2 

Substitution of (12) and (13) in (8) and (10) and equating powers of £ , 
the following sequence of equations are obtained 

4- b, (e,-T.) - 0 (14) 

at ' 

+ + b, ( a, -77 ) =<* (is) 

a a t * d 

+ + =: * • , i>Z (16) 

H f zrC- 

for the fluid and 
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25 



— bi. ( & a -To) - o 



(17) 



- b 2 ( 0,-77 ) = *, 

7>t ' 

for the solid. The boundary and initial conditions are 



(18) 



<->2. (19) 



and 



1 


° J 


*fc> D 




0 4*1 <ex> J 


t zs 0 


(*1j t) =-0 


0 J 


£> 0 ^ 


104 (»1^ ) = ° 


0 4 1 <co 




T 0 (*i^) 


1=0 j 


tz° 




0 <’1<« j 


t z= O 


=0 




t £ o 


TiM) = o 


°=i < °° > 


t = o 



( 20 ) 



,;>i (2i) 



( 22 ) 



6i (23) 



Since the order of these equations is less than that of the original 
equations, only one boundary condition can be imposed. It is not 
expected that the solutions will be uniformly valid. Starting with (14) 
and (17), we take the Laplace transformation with respect to t. The 
transformed equations are 
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+s + b, (e*-fo) 



— 0 



(24) 



( S+ bi )X = 



(25) 



Expressing T 0 in (25) in terms of Q g , substituting in (24), and solving 
the first order homogeneous differential equation yields 



T, (v*) = 



S(s+ bi) 






-1^ (^S-h bl - 



bi b?- \ 
S+ b 



(26) 



Ml > 5 ) 




-1 (s + b, - 



bi bz. *\ 
S + b 2 / 



(27) 



The procedure for inversion of (26) and (27) is as indicated in 
Example 7-5 of [2]. The results are 



Bo (l;6) =0 



for t<H (28) 



0* = e 



-V 



r 



• T » ) + b * I 0 [z {ay ) j 



for t> ^ (29) 



and 



X W’*) = ° 

t 

b z Je bz ^l 0 



0 



for t<1 (30) 

for *>*1 (31) 
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where a =■ anc * £ = t- \ . Note the use of a in this sense; it 

should not be confused with a linear dimension. It may easily be 
shown that Q 0 (*l,t) and To( *1 , t) both tend to unity as t -+-oo . 

It is clear that the two solutions obtained for T c ( q. , t) and Q 0 ( q, t) 
have a wave like character with no disturbance for \ ^ t. It is well 
to pause a moment to refer back to the physics of the problem. For 
fixed n. (i. e. , fixed value of x) no disturbance takes place until t > v^, 
or for fixed t, there is no disturbance for the region ^ > t. In what 
follows the condition T = t - H = 0 is referred to as the "wave front. M 
The pipe wall temperature T 0 ( Vl , t) is continuous at the wave front 
with discontinuous slope, while the fluid temperature Q 0 ( q , t) is 
itself discontinuous at the wave front. 

The first corrections 0, ( *1 , t) and Tj ( , t) are governed by 

equations similar to those for 0 O ( , t) and T 0 ( t\ , t) except that these 

new equations are nonhomogeneous with the second derivative of 
0 O ( *1 > t) and T 0 ( , t) as forcing functions . These non-homogeneous 

terms vanish ahead of the wave front, and thus no correction ahead of 
the wave front is expected in the sequence of solutions , T 0 , 0*, 

Tj, ... . In the next subsection we will demonstrate this by actually 

constructing the solutions 0^ and T^ and will see that the present 
perturbation expansion fails at the front. 
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B. HIGHER ORDER SOLUTIONS 



The corrections 0, ( ^ , t) and "J^ ( n. , t) are also obtained by- 
using Laplace transformation. From (18) 



v ' m H 2 - 



( 32 ) 



and from (26) we have 



aV. _ t. 



’brf S(S+b 2 ) 






(33) 



Likewise from (15) 



M + 0, 

bn. 



= biTi + 



* bT J 



(34) 



and from (27) 



iL. = i(s + b,-LL_f.&<pf-n(s + b,--bh-)] 
^ S V St bi > r [ L V S-rbi'J 



(35) 



Substituting (33) in (32) and solving for T^ in terms of Q t and then 
substituting the result in (34) results in a first order non-homogeneous 
differential equation for 0, (V^, s). 




(5 + b,/ 





-\(s + bi 



bt_biN 

S+bj' 



(36) 



The solution of (36) is easy to obtain. 
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(37) 



§ ' =5 - n lw + ^"i77 t ) ^ P H 



S-+- t»i - 



bi bi 

s+b 2 



)) 



It is now merely a question of substitution of Q ^ from (37) and 

yfo 

we have 



from (33) in (32) to obtain T^( , s). After rearranging 



*t~ S bi- 
ll m 



( s + bij 



^ J<m_ 

l* L (s+bi ) 1 S 4 bi 






(38) 



To invert (37) it is convenient to express 0^( ^ , s) as follows 



8 , = s e 



* a 4 

-n bl -ns rr^j, , . : n 

e • e • e ) (s + bj.) 

rt- 0 



(39) 



where 



k = <x n 

o f v 

kj = 2* b,1 

k 2 = o< m b|bi ^ 

k 3 = 

k 4 = ^ 



(40) 



Similarly 



a S 

= -lb, -VJ 5 s + fex 

T = s e e -e 



(W 



n=i 



(41) 
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where 




( 42 ) 



For the inversion it is useful to employ the formula 




(43) 



for integer n > 1. For n = 0, the formula is made valid by adding 
the term S (t) to the right side; this is the Dirac delta function. In 
order to express the result in terms of modified Bessel function of 
first kind order zero and one only, we use the following recursion 
formula 



for modified Bess el functions (Cf.W. G. Bickley [5], Number (20a)). 
It is also found convenient to introduce the function 



I n+1 (x) = I n _i(x) - (2n/x) I n (x) 



(44) 




(45) 



for which the following differention formulas are valid. 




(46) 



In order to invert (39) consider the following function 
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a 

s 



fM= e s ^\s" 



F(i>*hf {M 



* -o 



(47) 



using the formulas, 



/ 



,-1 



= S(t) + ff X-i [i ) 



(48) 



and 



/ 



a 

-n f 

S.e 






(49) 



yields , 



F(n.t; = -uw+>t. (^)"'r„-,(^) 



yi 



(50) 



Next using the recursion formula (44), (50) may be expressed as 



F( \ , t) = k Q £ (t) + m 0 I 0 (2 ^aT )+m^ (51) 



where 



f m o = kj + (Va) (k - ilL2_) 



lm 



, = *1 <* f ( a /t) + k 2 + 2 - _k3_ +k 4 (V 



(52) 



a) 



Thus we have 



l 



-I 



± ± 

S+bj 






n = o 



-kftr 

e Isw+Wol 0 (2jat) + 



(53) 
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Now let 



G (»\^s + bi) = s * f ( 1 ,5 4- bz j 



(54) 



Thus 



r\r 



— J <1 I * ^ 

G{\^) = £ a( r lj^b i j|=/ S-f(l,S+bi) 

= ^ F(^> + j| + / ' ^ f S + b 2 )- F( 0 + ) 

-4s(t) 4 - dp ft,*) 



To perform the indicated differentiation, (46) is used; thus 

V 



(55) 



^(vtj = e 



where 



oc^S(t) + ^^(2i>,-bi)S(t-j + C t I 0 (i{a+)+ (56) 

ql^ijat)] 



Co- ■ k - l ^yi^f ^b>^ f (b,.2bjt0C m b i + ^b,t (i>,-zb 2 + |) 
c i= - {% A«ba)(a/t j Vz + ^ i ( 2 b, - b L ) (a /t J Vz + 



(57) 



!4 



( 2 * lj b*1- * b* -tf b, - * b 2 ; b 2 ^ - 2 * b, fa + 2 * b, b,) (t/a ) Z '- 



(*m b^ Lz *1 J (f/a) 



Vz 



It is now easy to obtain the inversion of (39), which may be written 
as 
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( 58 ) 



-n b * -v _ . 

0, (l^> s ) = e • e • g (u s+ bz) 



Thus 



0,(^,t) = o 



-tX b| 



e,(u*) = e GCn^-i) 



for "t < 1 

for t> 1 



(59) 



proceeding in a similar way, consider the function 

JkX 

A„ . 



f*(i»s)= e s £C-s n 



F 



f (i'S) 



Thus 



F*( V t) =p: (±) 2 i,.,(2/5t) 






(60) 



(61) 



Expanding the sum in (61) and using the recursion formula (44) yields 

L(zffi) + f % (62) 



where 



¥ 

yn 0 



it 

vn> = 






(63) 



Again let 



G*(>\,s+b2) = S f *( r l 5 S+ b 0 



(64) 
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Thus 






fa 



= / 
-I 



-I 






= / +- /' Sf^^.svbij-Ftv®*)] < 65 > 



= 4* Stt) +i F*(vtJ 



The indicated differentiation in (65) is obtained using (46). Thus 



G (\,t) = e 



-bit 



* b 2 ^(t) + C I 0 {z{at) + cU i{z( rt) 



( 66 ) 



where 



c? 



(°< m b 2 +- 2b| bz *. t>2 1 J + (t/ a J 0<*b|b2 , 

. (2b,t> 2 ^4-^ b, - b 2 ^ + lb 2 + -b|bzH°<f/o< ) 

^ 3 i 

('t/a) o< m b, b 2 t. 



(67) 



y 

Cj — ( £J *+. 0< m (b,b\^- bi/rj- 6b,b i /r|-2b2/ K j - 8 /^ l + 2b^b 2 J 



+ * b,b t ^ ( b, -2 b 2 -*- 1/^) (t/a) 



Vl 



+ (°$, ^£ , 'i)(b,-zbi.-Vr l ) [t/a) 



?/z 



It is again easy to obtain the inversion of (41), which may be written as 



— -Y[S 

h (l>s) = £ . e • I (a ,s+b 2 J 



( 68 ) 
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Thus 



Ti(v t) = o 

x -bil * 

Ti (v *) - e • 6 (*l> t-»\) 



for 

(69) 

for t z r i 



Although the inversion of , s ) and T^( ^ , s ) were performed in 

as compact a manner as possible, it is desirable to rewrite 6 t ( , t) 

and T v\, , t) as 

=T7(vt) =<> for t <*l (70) 



-b,'| -bit \ / 

ftOv*) = e -e ^ f lS(tJ + (zb|-bi) S C^) + 



C olo{2-{crf) +C|Il(2(a?)J for t>ri (71) 



and 



-ri a 

Ti(\>9=e .e 



+ CI„(?/ar)+^T,(j/af)j 



(72) 



for t>*i 



It is of interest for future discussion to obtain the limiting values of 
9j( *| ,t) and Tj ( n,t) as IT — * o + 
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^ => £. £yrx 

T-*-O f r-*-o + 



-k.t * 

;r 6 

r~>o + 






CLU^ + c,!, (2 far; 



^ f rb|^ ^ (bi-ZbjJ+«wb2 



(Itar 



-o^a^ 1 c< f rt (ibi-bi)^ 

i>bzrH + fF 



ff-rfl...' 



Thus 



*&- 

tU<f 



ft(-l.t) = 



-b,t . L 
£ . o< bi 02 



\ bibit^ 2 ’+^ (zbi-bi) 



■ bl y . °(t» 

+ bi" 2 + «( 



And similarly 



(73) 



-b,7 



Ti Co I„ ('lie 1 -'? ) + C j If (z 

r-> o + r_ ^° 



-b,*l 



= e 



Z , l 



t* 2 \ + <* b 2 bl Ij 



(74) 



Since the solutions obtained in (71) and (72) are complicated, they are 
organized so that the character of the solution can be seen as clearly 
as possible. Equation (70) shows that there is no disturbance ahead 
of the wave front. It is also noted that at the wave front (T = t-7= 0) 
0^ and are singular, both containing delta functions and 0 j(*l ,t) 

being dominated by the derivative of the delta function. Thus 0i(>T , t) 

2 

is not of order unity and £ 0| is not small in comparison to 0 ( ^ t) 
at the wave front. 
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Although it is evident that the labor of obtaining higher order 
corrections (that is, higher order terms in the perturbation expansion) 
would be most formidable, it should be equally evident that such 
additional corrections would have a nature similar to that of the 
corrections already obtained; namely, they would exhibit discontinuous 
behavior at the wave front and would fail completely ahead of the wave 
front (that is, they would be consistent in predicting identically zero 
temperature). However, it is to be expected that they would possess 
calculable limiting values, from the left, as the wave front is 
approached. In effect, we thus have what appears to be a perfectly 
valid way of determining the behavior behind the wave front, but one 
which fails at, and ahead of the wave front. The nature of the failure 
of our original perturbation expansion is characteristic of singular 
perturbation problems [ 9 ], and we are therefore led to seek a singular 
perturbation expansion which will resolve the difficulties which have 
been encountered. This is done in the following subsection. 

C. SINGULAR PERTURBATION ANALYSIS 

We have seen that the perturbation solution represented by the 
series (12) and (13) involve no difficulties, except those of analytical 
complexity, for any > o , that is, at any point behind the wave 
front. However, this solution is useless for points ahead of the front 
f <0 , since it predicts zero temperature for fluid and solid at all 
points ahead of the wave front; physical considerations indicate that 
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diffusion through the wave front will result in non-zero temperatures 

ahead of the front. A mathematical indication of the difficulty may be 

2 

seen from the fact that £ is not small at the wave front, and 

4 

presumably the same thing could be shown for £ q ^ , etc. Accord- 

ingly, we now look for another solution which will be valid ahead of the 
wave front. To do this it is convenient to introduce the new variable 

Y = (*t -t) / 6 = x - Vt (75) 

and we now seek new solutions 0 ( ^ »t), T( | ,t) which will be valid 

for * 

Making the change of independent variable from ( ^ , t) to , t) 
provides the differential equations 



o< 



~b~B 



i Z^ 2 






fc.(e-T) 



(76) 



- a-; at 



£ t>» ( e - t ) 



(77) 



Note that the transformation (75) changes our viewpoint so that, in 
effect, we are now moving to the right with the front so that the fluid 
appears to be stationary and the pipe appears to be moving to the left 
with velocity V. 

Note that in contradistinction with the situation dealt with pre- 
viously, now g (to the first power) appears in the governing equations. 
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Accordingly £ , not £ , appears as the natural parameter with respect 
to which we hope to construct a perturbation expansion. Thus we deal 
with 



9 CM) = 0 O fiA) + £ 01 (Vt)+ f (78) 

= To [^A) + € Tt i^A) t (. Ti fa A) +•" (79) 

(Later when we "match 11 these new solutions to the solutions discussed 
in the previous subsection, we will add an asterisk to distinguish the 
new solutions. 

The conditions that all these new functions must satisfy are 

0^ v o j Til — ** 0 as ^ (80) 

0/ = 71* = o for t - o . (gl) 

for i = 0,1,2, ... 

Additionally, they must satisfy appropriate conditions at the wave 
front. Actual selection, or determinations, of these conditions for 
f = 0 will be made in the next subsection where the matching process 
will be treated. 

Substituting (78) and (79) into (76) and (77) comparing powers of 6 , 
gives rise to the following sequence of equations. For the fluid, we 
will have 
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(82) 



+ b, (e.-T.) 

e* = e, + b, (e, -HJ ( 83 ) 



* 0* =4 + b ' f & 2 -Tz) 

• • • 

and for the solid we obtain 

To 



(84) 



(85) 



' . 

T — ~ 1 1 4* l 0 - b Z (®o”7o ) 
T = -ii + 7T - b z ( ^ -r;; 



( 86 ) 



(87) 



(where the prime denotes partial derivative with respect to ^ and 
dot denotes partial derivative with respect to time). 

Condition (80) and equation (85) show that 

Tc= 0 (88) 

Now consider the first non-trivial equation in the sequence, 
namely (82). This is a classical diffusion equation, the solution of 
which, subject to (80), is 

-b ,t 

3>(T' f ) = A„e ■^fc(t/ 2 frt) (89) 
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The constant A Q is available to satisfy a boundary condition for Tjj =0. 
The next equation to be considered is (86) which in view of (88) may be 
written as 




Together with conditions (80) and (81), we may write down the solution 
of (90) as 

bt °° 

= A 0 o 2 e ' Jerfc ft/ 1 fig) ^ (9D 

X 

Note that no new constant of integration is required in the last step. 

The next equation is (83) which is similar to (82) except that T^, 
given by (91), is involved. The solution could be obtained but we will 
not do more than remark that another constant of integration will be 
introduced. Then (87) provides T£j no new constant is introduced. 
Next (84) provides a means of determining Q 2> another constant of 
integration is involved. 

Briefly the solution for Q ^ (i = 0, 1,2, ... ) involves a constant 
A^, which is actually a function of time t. Since the equations for T^ 
are first order, their solutions do not require the introduction of 
arbitrary constants or functions, equation (80) and (81) providing the 
necessary conditions. 

No difficulties are encountered (except those of increasing 
analytical complexity) in continuing this process and we expect 
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convergent series (78) and (79). However, as yet we have made no 

move to relate these expansions to any conditions at the left which 

might serve for the establishment of the as yet undetermined functions 

A (t). We proceed to do this in what follows immediately, 
i 

D. MATCHING OF SOLUTIONS 

-JU JL. 

T 'T' 

We now employ an asterisk to distinguish the solutions ( 0 and ) 
of the immediately preceding subsection, valid ahead of the front 
(i. e. , for ^ > o ), from the solutions ( & ^ and T\) of the earlier 
subsection, valid behind the front (i. e. , for 'C > 0 ). We will employ 
the following notations. 



c£i-m ( x , t ) =. ( t ) - di 

X-+Vi~ 



(92) 



V (*, * ) 'Tc (t) = Ti 

X-*W 

[^*(v)j - t *■ j =0.* 

(Tift.*)] ■ T. (t) =T. 



(93) 



(94) 



(95) 



(Note that here an overbar does not indicate a Laplace transform. ) 
Examples are 

0 O = e' b,t (96) 
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0 



( 97 ) 



0i 



T, = e 



e" b,t .o< U 2 tU l>,b 2 t Z + (2b r b 2 )t-+ k-Z+ 

t I t>z 

b,t (*>i + 2b|bz* * -*b\t + * b\b,t Z J 



= 0 



-b,t 

A a e 



CD 

= A , b z e^Jerfc (t/ z fiTi t ) d£ 



( 98 ) 

( 99 ) 

(100) 

( 101 ) 

( 102 ) 



o' J ; 

The functions 0^ and T . contain as yet undetermined functions of t, 
namely A q = A Q (t), A^ = Aj(t), etc, , and we must be able to evaluate 
these systematically. Our procedure is to evaluate so that the per- 
turbation sequences to the left and right of the front, truncated so as 
to include only terms containing € to an arbitrary even power (and 
lower powers) will match at the interface. 

The first matching is to require 



- - * 

6o = 00 (103) 

To =T* (104) 



Equation (104) is automatically satisfied. Equation (103) provides a 

condition which can be used to determine the constant A . We find 

o 
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A Q = A Q (t) = 1. Thus, our first (matched) approximation is given by 
( 29 ), (31), (89), and ( 91 ) with A q = 1 in the latter two. 

The second matching is 



- i - -* -* * - * 

+• € 9, — 9 0 + £ 6) + £ 



(105) 



— l - -* =. * 1=.* 

To + e r, = T 0 + + £ T 2 (io6) 

^ 5{< vu ❖ 

The functions 0i , 0£ , ^1 » an< ^ ^2 ^ ave introduced two additional 

unknown functions Aj(t) and A^(t). Equations (105) and (106) provide 
conditions for their determination. 

Similarly the next matching is 



1 - 4- 

0 o -h 6 0| + 6 



— ^ 2 — X 3 

0 O +■ € 0| + £ ©* -4- £ ■+ £ 0^ 



(107) 



T 0 \ tT, + 



H — 

* Tz 



— T* + £ f 6 Ti + £ T 3 -+ £ Tlf (108) 



in which two additional functions, A^(t) and A^(t) appear and which 
afford means for their determination. Proceeding in this manner, all 
functions A^(t) may be determined. At each stage we have matched 
solutions to the same degree in £ , on either side of the front. 

By this procedure, eventually (1) and (4) are satisfied. Moreover, 
at any stage, the necessary coefficients A^(t) are calculable so that 
the approximation affordstruncated series , available both to the left 
and to the right of the front, with matching at the front. Thus it appears 
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that we have succeeded in describing a process which provides two 
sets of series expansions, one valid behind the front and the other 
valid ahead of the front. (A ,r et" consists of a solution 0 for the 
fluid and a solution T for the solid. ) The labor of constructing the 
successive functions that occur in these series is considerable and 
we have actually obtained only a few of them. We have actually made 
the ,! first matching 1 ' and could, with additional labor, determine 
sufficient additional functions to accomplish second, and even third 
matchings. However, the writer confesses to a feeling of uncertainty 
in this entire matter of constructing the ordinary and the singular 
perturbation expansions and the matter of performing a matching 
where it is necessary to join them. Accordingly he suggests that it 
may be fruitful to consider this matter further, employing the insights 
which might be provided by a more sophisticated mathematical 
treatment. 

E. A DIFFERENT SOLUTION AHEAD OF TEMPERATURE FRONT 
The solution given in this subsection may be of some interest 

9 

i 

although it involves a certain basic inconsistency. In this solution 
we are concerned with the temperature variation ahead of the front 
and are able to consider both radial and axial variation in the fluid, 
as well as both axial and radial diffusion in the fluid. However, we 
assume that the temperature of the solid is uniformly zero and that 
there is perfect heat transfer through the interface; that is, we 
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consider an isothermal boundary condition for the fluid. (It would 
not alter the problem profoundly if a finite surface coefficient of heat 
transfer were assumed. ) Thus, this problem is not a conjugated 
problem. For a left end boundary condition we assume that the fluid 
temperature is a constant, independent of radius, which depends on 
the location of the front. Accordingly, the degree of approximation 
and inconsistency in this solution is considerably greater than is to 
be found in other problems treated herein. However, the results may 
be useful in considering the temperature distribution ahead of the 
front at quite small values of time, before significant radial variation 
has built up in the front. 

We will consider equation (II- 1) and (II-2) with isotropic thermal con- 
ductivities and subject them to a change of variable 

§ = x-Vt ; t~t (109) 



Then equation (II- 1) reduces to 



■ff I'M ,qe _ j_ be 
r br br 2 - ^ bt 

Similarly equation (II-2) reduces to 

^ 1 ** dr br 2 - -dt 



( 110 ) 



(ill) 



We note that the effect of the change of variable, (109), is to change 
our view point as follows; the fluid remains stationary and the pipe 
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wall moves to the left with the same velocity V of the fluid. If the 
velocity V is large then we state the following assumption 

Assume that for T£>0 the pipe wall is essentially at zero 
temperature. 

Hence we will concentrate our attention on equation (109) and 
find the solution when equation (109) is subjected to the following 
conditions 
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condition (114) yields 



o <r<a 




t>° 


(112) 


o < r < a 




t> 0 


(113) 


r = a 


/ 


t> 0 


(114) 


r =0 




t > o 


(115) 


o <r < a 




t = 0 


(116) 



— ► oo implies x — =► 00 ) 



«»\) (in) 






( 118 ) 



and condition (115) is already satisfied. 
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Next, taking Laplace transforms of (110) and (117), we find 




A solution of (119) which satisfies (113) is 



i 




Thus, writing , we have 



§ “ £c, 




Condition (112) then yields 



oo 

C J. (v) 

n=l 

Using the orthogonality of Bessel functions, viz. 

ft 2 

Jj. (^.r) J a ( X„f) «• Jr = Y J- ( x » a ) 

0 

if n = m and zero if n ^ m, we find, 




( 119 ) 



( 120 ) 



(121) 



( 122 ) 



(123) 



Hence 




(124) 



9 - z*}' ^ 

M\f>Ml*» a ) 

n— I 




(125) 
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Inverting (125) we find 




-\Sf 

e -€rfC ( 







(Note ; that solution (126) is for the case when the wave front is 
maintained at unity. ) It is now possible to use Duhamel's super- 
position integral when the condition at ^ = 0 is 

, P x - b ‘ t 

Q 0 [\>*) = e 



(127) 



which is the simplest estimate we have, Cf. (96). 



e(1,r,*) = 










. d$ 

2=t-S 



Thus, we obtain 

(128) 



where the term in brackets may be obtained by differentiating (126). 

Obviously, this solution is for cylindrical geometry, but the case 
of slab geometry differs only in detail. 
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V. CASE WHERE FLUID TEMPERATURE IS AVERAGED AND 
AXIAL CONDUCTION IN CONDUIT IS NEGLECTED 



A. SLAB GEOMETRY 

1. Development of Solution 

In this section equations (11-19) and (11-20), together with the 
appropriate boundary and initial conditions will be further developed. 
For convenience of reference these equations are rewritten below. 



V T>9 _ j_ 


* i 


(0-T) 


(1) 


dx l 'ftx <*, 
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ii 
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(2) 


0 = 0 


X < «> 


, t =» 


(3) 


0=1 


X = o 


t t > 0 


(4) 


0-*- 0 


X— ► 03 




(5) 


T-° 


0 <3 < b 


, t = o 


(6) 


II 

O 




^ t 


(7) 


*3 






-k£M ( 9-T) 


3 = ° 


j t>5 


(8) 



Note, that 9=9 (x, t) appearing here is the average fluid temperature 
as defined in (11-18); we have omitted the asterisk uniformly throughout 
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this section. Also, T = T (x, y, t) is the temperature in the duct wall; 
we have omitted the subscript 2 on y; in this section, the symbol y 
denotes y ^ in figure lb. We also note that this analysis neglects 
the axial conduction in the duct wall. Finally, note that the surface 
heat transfer coefficient appearing above is the ’'equivalent 1 ' film 
resistance coefficient £ discussed following (11-17); we omit writing 
the tilde in (1) and (8) above and in what follows. 

In order to solve (1) and (2), Laplace transformation is used. 
Equations (1), (3), (4) and (5) become 



to V ^ 

lx z ~ ax 



s 9 + — 



oi 



f 



a K 









X = 0 



Q 



00 



(9) 



( 10 ) 

( 11 ) 



Similarly equations (2), (6), (7) and (8) become 



tf = if 






o <x <.oo , u = <? 



> 




( 12 ) 



(13) 



2T=o 



O ^oo 



(14) 



The solution of (12) is 



w 



T = A e 



B? 



(15) 



where 




(16) 



Condition (14) permits evaluating B in terms of A, and we arrive 



at 



T = A 



fa $(^-3)) 
e + e 



(17) 



Since the terms in parenthesis are clearly independent of x, the 
imposition of condition (13) implies that A in (17) must be a function 
of x and s. Thus using (13) to express 0 in terms of T, we find 



9 = (i (5) • A (x , s) 



(18) 



where 

Substituting from (17) and (18) into (9) yields a homogeneous 
second order ordinary differential equation for A(x, s). After re- 
arranging and substitution for (s) this becomes 

- Y M = <](sj. A (\jS) ( 

bx 1 d 
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where 



S_ + j_ / K.{W( H) \ 

“V a K f \ I+- K»; 

The general solution of (20) is 




A(XjS) = C-e* p 





( 22 ) 



where C and D are constants of integration which depend on S. The 
boundary conditions (10) and (11) in conjunction with (18) define the 
following conditions for A (x , s) 

{x ' s)= x= ° <«> 
A (x,s) >0 X — > OO (24) 



Condition (24) clearly implies that D 5 0 whereas (25) yields 






Therefore 



A(x,s) = 




(25) 



(26) 



Substitution of A(x, s) from (26) above into (18) and (17) yields the 

fluid and pipe wall temperature distribution in Laplace transform form; 

viz. 
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(27) 








i fipfaJUm V 



and 



T(x,^s) 



✓ \ 



c*sl[%(b-y)] I 

c °^ (b§) + Jip(b pA^&ibpj 







(28) 

f(^)-0(x,s) 



In order to invert (27) and (28), we first consider 0 (x, s). Denoting 
the radicand by m we have. 



m = 



( y./ + i, W w (b 1 )[K ? (X) 



Kr 



• ■ 3- 4 

bV * t - abKi J 



Km 



l+ ^(b j)W(lf) 



(29) 



It will be convenient to cast this equation in another form. Write 




so that 

s = -« X* ; y a 

Then (29) becomes 

_ bX ta*. (b\) 



(30) 



(31) 



(32) 
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where 




j 




(33) 



are Biot numbers for the wall and for the fluid, respectively. We 
are concerned with finding the singularities of (27). Except for the 
pole at s=0, these will all be branch points corresponding to the zeros 
and poles of m . Considering the numerator and denominator of (32) 
separately, we first note that the denominator has poles for bX = 
(n+V2)tf , n = 0, 1, 2. . . . These will also be poles for the 

numerator unless 



We may assume that this equality does not hold for any integer n 
for otherwise we could make an infinitesimally small change in any 
of the physical parameters so as to assure that the equality is violated, 
without making any significant change in the physical problem. Thus, 
the only zeros and poles of m correspond to the zeros of the num- 
erator and denominator, respectively. 



The condition for zeros of the numerator (i. e. , zeros of m ), is 

that 




(34) 




(35) 



where 
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( 36 ) 



/ V \} 

B* \ Zo< f / " A 




and the condition for zeros of the denominator (i. e. , poles of m ), 
is that 

tan(b X ) = f 2 ( X ) (37) 

where 

f 2 ( A ) = JLn»- (38) 

bX 

The real zeros of (35) and (37) may be found easily. Their location 
is illustrated in Figure 3. The relative position of the dashed line 
and the dot-dash line in Figure 3, is important. For small X it is 
evident that f^( X ) ^ ^ )• Also f^( X ) has a vertical 

asymptote for a value of X greater than that of its single zero. Thus 
the positions are as shown in Figure 3. 

Finally, we must assume that there are no other zeros than the 
real zeros indicated in Figure 3. However, both (36) and (38) satisfy 
the hypotheses of Appendix B which shows that there are no non- real 
zeros . 

Thus it is clear that the only singularities of m are poles; 
furthermore, m has only real zeros and poles, infinite in number, 
and alternating with each other, of which that having the smallest 
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y = tan (bX) 

y = fjlX) 

y = f 2 CX) 

o - zeros, x - poles 



Figure 3: Poles and Zeros of m(s) 

(Slab Geometry) 
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numerical value is a zero. If we designate the abcissae of the inter- 



sections denoted by 0 in Figure 3 by the symbol \ and those 



n 



denoted by x by the symbol , we conclude that the singularities 

of (27) are: 

a. branch points (corresponding to zeros of m ) at 



= - o< \ * 

^m A n ’ 



and 



b. branch points (corresponding to poles of m ) at 

s = - cJ *f 

m / n 

c. a pole at S=0 with residue 1 

Furthermore, the branch points of type a and type b alternate with 
each other, that closest to the origin being of type a. 

Thus the hypotheses of the inversion theorem given in Appendix A 
hereof are satisfied. Accordingly, we may immediately write the 
inverse as follows 



0(x, t) 

where /Ma) = (| m ) 



b K 

CO f -'l*' 

= 1 ' J \ ^ 

K=l \ 

, . . . , 



<h 



(39) 



Similarly T(x, y, s) as given by (28) is conveniently organized in the 
form of the product of two Laplace transforms. Accordingly, its 
inverse may be easily obtained in the form of a convolution integral. 
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( 40 ) 



CO 



T(x M +) = j Y' — — — ^ Ant> ^ Co% 

9 L 2 > n ^/^(zx n b) 

M=f _t 






\ t 

mA n ^ 

e 



0< \ 2 1 
m"n l 



£ ” " • 0(x,^)d»X 



O 

where X* are the roots of (37) and the integrand contains the 
function (39). We have used the formula 




Cos£ 

Cosi (bs'4j + i ( b&) ( bsty 

4 



4 



\ X M ^ " A^-rt (XnOlJ CoS (Anb) 

Z--1 2 X„ b + ACr\ (2bXnJ 



\ 



(41) 



which may be derived directly or inferred from equations given in 
[2], examples 5-3 and 7-11. 

2. Asymptotic Approximations 

Thus the desired solution of the system described by equation 
(1) through (8) has been achieved. We will not in this work discuss 
computational aspects, but will now obtain some asymptotic 
approximations . 

Equation (27) can also be written 



e(x,s)^ I e*p 



T <- * ^ ~) 



(42) 
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where 



X*I^W(bj) 



ol 2 C Z L 4 5,6 4 .S 

Si _ s o 2 . s k _ in s' b 

°<m 3«* + 7? ^ 35 ** + “ * 



(43) 



Consider the radicand to be expanded in ascending powers of s. 



e(x,s) = -Lc* P 





(44) 



Considering small values of s, which corresponds to large t, 
the first asymptotic approximation is obtained by neglecting all powers 
of s except the zeroth power. Correspondingly 

§(x, S) = i (45) 

which yields 

e (x,t; = i (46) 



This is evidently correct as an asymptotic expansion since for any 
finite x the temperature approaches unity for sufficiently large time. 

A nontrivial asymptotic expansion is obtained by retaining 
the first power of s. Correspondingly 






£-x, 


/fy. ) l + 
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/A>vf 
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4 



(47) 



The reason for writing A m /Af , the area ratio, in place of 
b/a which naturally would appear in the last term of the radicand will 
appear in Section VB2. 
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The inversion of (47) is tabulated [29], thus 




( 48 ) 



where 






(49) 



is a retarded velocity which is identical to that obtained by Munk [24], 
and 



is an effective thermal diffusivity for the fluid. This solution should 
be compared with that given in (VII- 10) where a different set of 
idealizations leads to a similar result. There, the retarded velocity 
(VII-5) is the same as (49), but the diffusivity (VII-6) is slightly 
different from that given in (50). 

For the corresponding solution for the conduit wall temp- 
erature we employ (48) in (40). 

Retaining the second power of s would yield a useful 
improved asymptotic approximation if one could find a simple explicit 
inversion formula. Unfortunately nothing simpler than the general 
inversion given in Appendix A has come to light, and this has provided 
the complete solution given in (39) and (40). 




(50 
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B. CYLINDRICAL GEOMETRY 
1. Development of Solution 

A similar treatment to that given in Section VA is employed. 
Only the significantly different details will be treated explicitly. 
Equation (2) will be replaced by 

tfr . Iff = ± ST ( 51 ) 

5>r‘ rX <,ji 

with the corresponding conditions 



T = o ' a <r <b , 


i = o 


(52) 


SI — 0 If* — fc j 

dr 


t > o 


(53) 


- K (o-t) r=a ; 

m <>r v 7 


t> 0 


(54) 



replacing (6), (7), and (8). Remarks similar to those following (8) 
are applicable here. 

As previously, Laplace transformations are made with 
respect to t. The solution of the transformed version of (51) is 

T = Al e ($r) +BKo{y) (55) 

where 

| (56) 

is the same as previously. Here I and K denote modified Bessel 
r 1 o o 

functions . 
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Proceeding as in VA1, one arrives at equations similar to 
(27) and (28) except that the function tanj} (qb) appearing therein is 
to be replaced by another function which we denote 



banB( - 



(57) 



(The initial b denotes "Bessel. M The remainder (anh) is intended to 
indicate a similarity to the function tan£(x) . ) 

Corresponding to (29), we have in this case, the equation. 

„ _ (58) 

_ \ (\A . d/L ) 



1 + J<a(bj) \,a~4 (if, <ViJ 



Corresponding to the well known relation 



4. tosn'd (x) — 



(59) 



it is not difficult to show that 



l b owi(Xj ft) — bcui ^X ; ft) 



( 60 ) 



where 



JW YAP*) -Wj'lpx) 



(61) 
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Thus corresponding to (32), upon defining 



X = _ C 



( 62 ) 



we obtain 




(63) 



Bn, - (bx) bom(bA,a/b) 



where 





(64) 



(Note the factor 2 in the present definition of B^‘ ) 



From this point, the development is precisely the same as in VA1. 
Considering the counterparts of equations (35) and (36), precisely the 
same treatment is applicable. Figure 4, shows a graph similar to 
Figure 3. Figure 4 was actually constructed by computer analysis in 
order to show that the function ban ( X , ) behaves like tan(x) ; the 

value of (3 used in constructing Figure 4 was fl = 0.9. . 



In the previous development, the fact that m had no non-real zeros 
or poles was proved by the argument contained in Appendix B. To 
construct a similar proof for the case in which tan x is replaced by 
ban (x, (i ) appears terribly formidable and we simply omit doing so. 

Thus, a solution corresponding to (39) and expressed in 
exactly the same notation as (39) (that is, identical to (39)) is 



There is one difference which we would prefer to ignore. 
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Figure 4: Poles and Zeros of m(s) 

(Cylindrical Geometry) 
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applicable in the present case. However, the function F (r,s) in (28) 
must now be replaced by 



F(<vO = 



LqWK.WtK^rjr.^b) 

(- f Wl a) + j K,(^)+ Kj'jojjl.qi .)' 65 1 



The inversion of (65) may be obtained using a solution discussed 
in Section 13. 4 of [7], by first specializing to the case of a unit step 
input, and then by differentiating the result with respect to time. The 
result, corresponding to (40), is 



OQ 




nn 




( 66 ) 



where jJ n are the roots of 

K *,/' 2 , » 



(67) 



and 



Jk(i ,x ) s (M) -Yk(( ' 68 > 



2 . Asymptotic Approximations 

It is possible to proceed here in exactly the same way as in 
Section VA2. The principal difference is that in the present case of 
cylindrical geometry, in place of (43) we now define ^ differently 
and arrive at different coefficients in the expansion. Thus 
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(69) 



X-(b|)b<u,e(l|,cr/b) = i (t '-<*)£ + 

where . . . denotes terms in s , etc. It requires painstaking work 
to establish that the coefficient of s° vanishes and that the first 
power term is as in (69) but the calculations are not illuminating and 
are not given here. 

As in the case for slab geometry, the first non-trivial 
asymptotic expansion is obtained by retaining the first power of s 
in (69); correspondingly, we arrive at equation (41) the inverse of 
which is given in (48). In writing (47) we anticipated the present 
development by writing A^ /A £ in the final term of the radicand so as 
to be able to accommodate to both slab geometry and cylindrical 
geometry with one equation. 
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VI. A GENERALIZATION OF A PROBLEM CONSIDERED 

BY ARPACI 



This case is described by (II- 10) and (11-14) subject to conditions 
that will be stated below. The problem description is thus 



_ V *© „ 
"bx* 7>x 
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(6) 



Q = Q (x, t), T = T(x, t) are average temperatures in the sense of 
(II-7), (II-8), but we omit using the asterisk which distinguished this 
significance. The heat transfer coefficient h in (1) and (2) above is 
defined by (II- 9 ) . 

This problem is a generalization of that considered in Section 
7-5 of [2], the generalization coming about by the inclusion here of 
axial conductivity in the fluid, a mechanism not considered in Arpaci's 
analysis. This problem is a specialization of the problem considered 
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in Section IV hereof since here we neglect axial conduction in the 



duct wall, whereas this mechanism is included in the analysis given 
in Section IV. 

Taking Laplace transformation with respect to time t, one obtains 

. d+b.ja + b,T = o (7) 

bx*- ** bx v *f 7 v 1 

ST = bz( e -T ) (8) 

where 

l -SP , -Kp 

b i = : b i - (9) 

' f>A„ 

Combining (7) and (8), one obtains 



ti 

1 x z 



V _ ( s. 
\ lx ' 



+ b,- 



bi bz 
S+ b i 




= 0 



( 10 ) 



subject to the conditions 



6 * J X ~ 0 (11) 

&— >0 X— >oo (12) 



Thus, the solution of (10) may be written 




(13) 



an increasing exponential having been rejected in order to satisfy 
(12). We write (13) as 
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(14) 



e = 



s r 



s 








2* f » 
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(S-f- bz ) 

4 



where 



A = 




(15) 



(16) 




LO - \>i + b) c<. 4 ¥—■ 
t 4 

It is easy to show that A>p> 0 • However, 

the magnitudes of X and p relative to b 2 . 
the function 



(17) 

we want to indicate 
To see this, consider 



f(s) 








(18) 



the vanishing of which corresponds to the roots p ^ X ; viz. 
f(- X ) = f(- p ) = 0. We note that the coefficient of the quadratic 
term is positive, and we note that 

f(-b 2 ) = -b!b 2 < 0. (19) 

Thus, we easily see that, 

A > b i > O (20) 

Thus the condition of the theorem stated in Appendix A are satisfied, 
and we can write the inverse of (14) as 
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( 21 ) 




Using the convolution theorem to invert (8), we obtain 

t 

T(*,-t) = h>z J (23) 

o 

We should remark that equations (1) and (2) above may be 
obtained by specializing equations (61) and (62) of Hiep's treatment 
[13]. However, we note again that Hiep’s results depend upon 
numerical inversion of Laplace transformations, a process which 
our experience indicates is of doubtful accuracy and reliability. 
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VII. AN IMPROVEMENT ON A RESULT BY W. MUNK 



Under the assumption that radial conduction is so efficient that 
essentially the same temperature prevails at all points of a particular 
fluid and metal cross section, W. Munk [24] obtained a very simple 
approximate solution for what he called "The Delayed Hot Water 
Problem, " In the present section we obtain a very simple improve- 
ment over Munk's result. Hiep [ 1 3 ] has also obtained the improvement 
given here and he mentions that it has been obtained by others. Our 
reason for presenting it is that it is a natural deduction to be drawn 
from developments already made in this thesis. 

We repeat equation (11-13) and the corresponding conditions. 
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We have used T to denote the common value of average temperature. 
By defining 



* 

v = v 









(5) 



( 6 ) 
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we may write 



o(*tL _v*£T = 

dx 2 lx. -frt 

Taking Laplace transforms with respect to time, we obtain 



( 7 ) 



**gf-v* 2 f = ST 
lx 1 £x 



( 8 ) 



The solution of (8) satisfying all conditions is 



T= - 



x /fv*Y\ ^ 

2** X VV2^7 



/* \2 



+ ** 



(9) 



Inverting (9) by the use of standard tables [29] we obtain 



T( x >t) = ' 



(1^) + e 



* Vj< * 

2 o<* 



er<c (*+!*) 



(10) 



This is the generalization described in the title of this section. 
Munk’s result is 

T (x, t) = H (x-yt) (ID 

where H ( ) denotes the Heaviside unit step function. Thus 

the present result (10) is like Munk's in predicting a temperature 

■v 

wave front moving with a reduced velocity V 9 but, unlike Munk's 
it also provides for a softening of the wave front due to axial 
diffusion in fluid and pipe. 
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VIII. DISCUSSION 



A. SCOPE OF ANALYSIS 

The general problem as formulated by equations (II- 1) and ( II- 2 ) 
for the case of cylindrical geometry and equations ( II- 3 ) and (II-4) 
for the case of slab geometry, even with the simplifying assumptions 
of subsection IIA2, are mathematically very involved. These 
equations contain two dependent variables, namely the fluid and the 
pipe wall temperature each of which is a function of three independent 
variables. There is also a formidable number of boundary and initial 
conditions to be satisfied. In Section III the problem is reduced to 
the solution of an integral equation the only unknown in which is the 
temperature at the interface, which is a function of only two variables 
no auxiliary conditions are required. However, as may clearly be 
observed by the reader, the integral equation of interest is math- 
ematically very involved, and retains a great share of the complexity 
of the original problem, for the following reasons. 

Firstly, the unknown function appears only under the integral 
sign (the so called integral equation of the first type). Secondly, 
there are two integral signs under which the function of interest 
appears. Thirdly, one of the integrals has a semi-infinite range, 
namely o to oo Fourthly, the kernels of integrations are 

complicated functions. 
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The complexity of the general problem, either in its original 
form or as compressed into this integral equation is such as to have 
led us to focus attention upon further simplifications which are more 
tractable from a mathematical viewpoint and the analyses of which 
are treated in other sections of this thesis. However, the integral 
equation formulation is believed to be of significance and its treat- 
ment, probably by numerical procedure, is considered feasible. 

Further simplifying assumptions are introduced in IID, in order 
to reduce the mathematical difficulties. In a search for the various 
methods which would render the general problem amenable to math- 
ematical analysis, it is found that by defining a suitable mean temp- 
erature in the direction perpendicular to the direction of flow, and 
by making a further assumption of constant surface heat transfer 
coefficient, the useful formulation given in (II- 10) and (II- 11) is 
obtained. 

Assumption of large velocity justifies the use of a perturbation 
technique, which is presented in Section IV hereof. 

The zeroth order terms, which corresponds to a solution given 
in Section 15-3 of [7] and also in example 7. 5 of [2],are obtained, and 
first corrections to the fluid and pipe temperatures are obtained 
explicitly. These solutions are found not to be valid at and beyond 
the temperature front. Accordingly a singular perturbation tech- 
nique is used to obtain an expression which traces the behavior ahead 
of the front. 
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In Section VI a further physical assumption, namely neglecting 
axial conduction in the solid, leads to a problem for which we are 
able to find the closed form solution by using the method of Laplace 
transformation. This affords an improvement over the solutions given 
in [7] and [2], cited in the preceding paragraph. 

In Section V the physical assumptions that are introduced consist 
of considering an average fluid temperature and neglecting axial con- 
duction in the solid. This case is solved using the method of Laplace 
transformation. Even though the methods of transform calculus 
reduce the partial differential equations to a system which can be 
readily solved in the transform domain, sometimes delicate math- 
ematical handling is required to invert the solutions back to the real 
variable domain. Such a case is that of equation (V-27) and (V-28), 
where the solutions in transformed form contain an infinite number 
of singularities (branh points). An elegant theorem for this inversion, 
established by Professor J. E. Brock, is presented in Appendix A. 
Thus closed form solutions are obtained. Furthermore, use is made 
of the properties of Laplace transform to yield some useful asymptotic 
solutions. 

In Section VII the view taken is that of W. Munk [24]. It was found 
possible to improve upon his solution by accounting for the axial 
conduction in the fluid as well as the solid wall. This treatment 
assumes, as Munk did, an instantaneous common constant uniform 
temperature at any cross section of the fluid and wall. Although 
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this result is obtained independently, it turns out to be known; see 
Appendix II of [ 1 3 ] . This solution is useful for small fluid velocities; 
Cf*Hiep [13], p. 32 and Preston [27]. 

We will make no further attempt to relate the problems considered 
herein to practical engineering situations. Such references as [7] 
and [2] are devoted to problems having similar relations to practical 
engineering situations. Their importance lies in the fact that they do 
provide definite mathematical solutions to well posed mathematical 
problems which have relation to the physical world of heat transfer 
engineering; it is the obligation of the person who employs such 
results to investigate the degree to which the mathematical conditions 
are satisfied in his own physical application. 

B . RECOMMENDATIONS 

Most important is treatment of the integral equation developed 
in Section III. Remarks about this are included in subsection IIIF. 

It would be useful to develop asymptotic formulas to employ in 
place of the exact inversion in cases where there are severe numerical 
difficulties associated with the latter. The method of steepest descents 
(e. g. Cf. Carrier Crook and Pearson - Section 6-3 [8]), appears 
promising for this purpose. 

There are several mathematical details which should be filled 
in (although in treatments similar to our own they are usually glossed 
over). These all relate to locations of the zeros and poles of 
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transcendental functions. In the case of equations (V- 5 9) and (III- 155), 
we have simply assumed that the zeros and poles are all real; under 
this assumption they are not particularly difficult to locate. How- 
ever a proof that there are no complex zeros or poles is really 
required. 

Inasmuch as instantaneously higher rates of heat transfer are to 
be expected at and near a rapidly moving front, it might be expected 
that an analysis employing finite speed of heat propagation would be 
enlightening. Accordingly, it is suggested that such an analysis be 
attempted in one of the simpler subcases dealt with herein; Cf. , for 
examples, Section 1-6 of [ 1 0 ] . 

Also since numerical difficulties are to be expected in applying 
the inversion theorem of Appendix A to the heat transfer problems 
herein, it is suggested that actual numerical calculations be under- 
taken. A multiple-precision facility, such as described in [ 1 4] and 
[28] seems necessary for this purpose. 

Another thing which comes to mind is to relax the specialization 
herein because of the assumption of slug flow, i. e. , constant fluid 
velocity. It seems clear that the assumption of velocity as a given 
function of the radial or transverse dimension would provide definite 
differential systems for analysis, but that these would, in general, 
be much more difficult to deal with than those considered herein. 

The treatment given by Hiep [13] relies heavily on the accuracy 
and validity of a numerical inversion of Laplace transform due to 
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Salzer [30], [31], [32]. However our experimentation with such methods 
of inversion showed that reliably accurate results could be obtained 
only in limited ranges and for relatively simple functions. This 
persuaded us to look for other analytical methods of treating the 
problem, such as perturbation method treated in Section IV hereof. 
However, it is probable that further development of techniques of 
numerical inversion may provide a useful and powerful tool for 
analyses such as those studied herein. 

We should also mention our feeling that no great difficulty would 
be associated with inferring response to a unit step change in left 
end temperature but retaining the other conditions specified in the 
analysis of Sparrow and DeFarias [33], who employed a sinusoidally 
varying left end boundary condition. 

Also although throughout we have assumed that the conduit had 
perfect exterior insulation so that there was no heat flux through this 
surface, it seems clear that other exterior boundary conditions (such 
as an isothermal condition) would, in most applications considered 
herein, involve only changes in detail and not essential changes in 
procedure. 
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APPENDIX A 



INVERSION OF SOME LAPLACE TRANSFORMS 



The theorem stated and proved in this appendix provides the 
Laplace inversion of (V_27), (V-28), and (VI- 14) of the text. This 
appendix is largely a paraphrase of [6]. 

THEOREM: If the function to be inverted is of the form 



f (*) - y 



-iMs) 



(i) 



V. J 

where m(s) is a meromorphic function whose zeros s = -a^ and whose 
poles s = -b^ satisfy 



Then, for 



0 < <C b| < - • • < ^ • 

t > 0, the inverse Laplace transform 



• < oo (2) 

F(t)= /V)] 



is given by 

F(+) 

where 



A r-V 



(3) 



Hi) =|/| m (-i)l' 



(4) 



The interpretation of the sum in (3) is as follows. If there is a 
finite number N of pairs a^ ^ b^., no special interpretation is 
necessary. If there is a finite number (N-l) of pairs a^ ^ b^_ and an 
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additional zero s = - a^, a N y t>N_ \> an pole s = - t> N = - oo 
is to be added so as to make N pairs. If there is an infinite number 
of pairs a^ b^» take N = c° 

In proving the theorem stated above we consider a sequence | R n 
of positive numbers, R n i a , bj^. for any n and any k, R n — >oo as 
n — boo, and we consider a corresponding sequence of closed contours 

C in the complex s- plane each of which encloses the pole of f(s) at 

n 

s = 0 but excludes all of the branch points at s = -a^ and s = - b^.. 

The contour C n consists of four subcontours C n i, C n 2, C n 3 and C n ^ 
as follows (see Figure A-l). C n i is a vertical line from s = - iR R , 
to s = iR n , indented as shown by a small semicircle passing to the 

right of s = 0. C n 2 is the arc s = R n 6 > H ^ Q <1T and C n4 is the 

<■ 9 

similar arc s = R n ^ • Finally C n 3 is the contour which 

passes to the right froms = -R n , along the negative real s-axis to 
s = - a j with small semicircular indentations into the upper half s - 
plane about each zero (-a^) and pole (-b^.) and then passes to the left 
from s = - a^ to s = - R n with small semicircular indentations into the 
lower half s - plane, as indicated in Figure A-l. 
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Clearly by Cauchy's theorem 



I„ = =e *p 



- \JnUo)' 



- exp 



-w 



(5) 



for n = 1, 2, 3 . . ♦ . Also, by the complex inversion theorem for 

Laplace transforms 

i..- <6) 

C*l 



as n — > 00 • 

To proceed further we wish to establish an inequality for arg 
[m(s)]. To do this we consider complex vectors emanating from the 
points s = - a^. and s ='b^, and we write 

u£ = 

&K 



Ytt 



**?{*■% I<) J 


; 


-ir<^ <tr 


(7) 


^pO-Yk) ; 


; 


-ir<Y£ <rr 


(8) 










W IT 0*£, 


) -r TT ( 




(9) 



= (S) (i/3n) 



(10) 



where 



(«> 



= pH TT u k -rTT* v k 



(id 



and 



(”) (n ) 

P* = 



( 12 ) 



The sums and products written with a superscript (n) indicate that 
only those terms are to be considered for which a n < R n or b < R n , 
as may be appropriate. 
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Now note for any s in the (closed) upper half plane we have 



t« S Yk > 



K = I jZ, 



( 13 ) 



N-> 

ft- -- ^TK-Yk) + %> 



(14) 



A geometrical argument is given in Figure A-2. Now either 
there are N pairs a^ , b^. each less than R , so that 

P. = £ ( *k" Yk) £ 0 

or there is an extra, unpaired zero, so that 

A/-/ 

_ 0 (15) 

K =: | 

In the first case we can write 



V-l 



ft- • *■ - - % J - Y„ 

i 4>, - Y„ < 4; < 'n' 



(16) 




Figure A_2: Demonstration that > -yj. 
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whereas in the second case we can write 



/V 

/i„= V JJYk-U.) £%<« 



( 17 ) 



K*l 



Thus in either case 

° < p* < fC (18) 

Similarly by considering an s in the lower half plane we find 

( 19 ) 

so that in all cases 

-rf< <rr ( 20 ) 

Thus, if we write 



i«« 



Si 









we note that 



< 21 ) 






(22) 



Rr 



- C °s(^"/z)] < -L 

Yt L J rvw 



uniformly on C n 2 > by virtue of the inequality (20). Thus, by Jordan’s 
Lemma I n 2 — >0 as n — voo. Similarly the integral on C n 4 vanishes 



as n — >oo - 

Next we consider the integrals around semicircles surrounding 
poles and zeros on contours C n3. The result is obvious near a zero 

si -ol k i 

s = • a k> e&e , U(s;|^ J- since > and 

I 1 a * 
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the length of the arc goes to zero as its radius does. Thus the integral 

on such semicircles about a zero, s = - a. vanishes with the radius. 

k 

However, for the integral on a semicircle about a pole, s = - b^, we 
must invoke (20) again. We have 



"l/m (s jj = e*p [- |A„ Cds ( < 1 



(23) 



even though — > oo as the radius of the indentation goes to zero. 

Thus the integral on such a semicircle about a pole, s = - b^, vanishes 
with the radius. 

We remark that although we have reason to believe, on physical 
grounds, the correctness of an evaluation by Kantola [ 1 8] , we feel 
that his argument intended to show vanishing of such an integral on a 
semicircle surrounding a branch point is deficient. Our argument, 
based on (20), is not applicable to his case. 

If now we write (5) in the form 



I„1 = I n - I„2 - I„3 - I„4 < 24 > 

and let n > oo , and also let the radii of all semicircular indenta. 

tions on C n 3 go to zero, we find 

F(t) = exp [- fi(0) ] - I 3 (25) 



where 




I 

ircc 




■e*p (si - y / rtUs7 ) <^S 



(26) 
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and now is a contour extending from s = - 00 along the upper edge 
of a branch cut, to the right, passing around s = - a ^ , and then 

extending along the lower edge of the branch cut to s = - oo . The 
branch cut is the negative real s - axis from - oo to - a^. There are 
no longer any semicircular identations in C^. 

To reduce (26) to a more convenient form, we first introduce the 
function (x), for real x, which denotes the excess of the number 
of zeros (s = - a^) over the number of poles (s = - b^) to the right of 
the point s = x. It is clear that 1^=0 in the interval -a^+i s 4 -b^ 
while V = 1 in the interval "b^ < s <£ -a^. Writing ^ now, rather 
than , since we have taken n — > oo and all zeros and poles 

are included in our accounting, we note that on the upper edge of the 
cut /3 = PfT/2 while on the lower edge ^ = — 2^/2. We also now 
write rather than jJ n . Thus 



x = -L 



1 e*p 



s-t - fj (i. 



As + 



'-00 

J-at 

-a, 



= l.l±e st 

ITU s c 



T 

£*p 



S -t - jj p (- i iZTf/ij d S , 



(27) 



-p t/ 2 j - eyp -p e*p 






-oo 



If V = 0, the quantity in braces vanishes, while if l) - 1 it 
becomes 

e* P (-ipJ _ exp(/ip) =-zi.AjZ+ix (28) 
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Thus 




( 29 ) 



If we now change the variable of integration from s to *i = - s, we 
obtain 



The use of the upper limit of summation is consistent with the remarks 
following the statement of the theorem, the proof of which is now 
complete. 

In equation (VI- 14) the solution obtained in Laplace transform 
form is a special case in which we have 




(30) 



and, thus, finally 




(31) 




(32) 



where lc > 0 < and o < Cl, < bi < 0Lj_ . The numbers aj, and 

bj are known explicitly by (VI-9, 15, 16). We have 




( 33 ) 
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and the inversion is 



-P© 

F(*J = e - 



l, 




r-it 




n A - 


p<*i) 


k 






( 34 ) 
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APPENDIX B 



NONEXISTENCE OF COMPLEX ROOTS OF EQUATION (V-36) 

Demonstration that such equations as (V-36) and (V-38) have only- 
real roots are occasionally called for in order to employ Cauchy's 
Theorem (and Jordan's Lemma) to evaluate the inverse Laplace 
transforms encountered in heat conduction problems. Carslaw and 
Jaeger [7] have some general remarks on this matter (Section 12. 8 
IV , page 324), making use of Sturmian theory; their reference is to 
Chapter X of [15]. Polya and Szego [26] consider similar problems, 
most of them concerned with the zeros of polynomials, but some.e. g. , 
Auf.172, of nature similar to ours. They use various methods of 
proof. Professor J, E. Brock has constructed the following elementary 
proof for the case encountered in equation (V-36). The case encountered 
in equation (V-61) wherein the simple trignometric tangent function is 
replaced by a generalization involving quotients of expressions involv- 
ing J- and Y- Bessel functions, is obviously much more difficult and 
we simply make no attempt to construct a proof in this case, taking 
refuge in the (weak) argument that the function we encounter is 
"similar M to the tangent function. 

Thus, we consider 

. c(«- s-* 1 ) ... ... 

xtanx = — : C>o , n>o o<b<» (1) 

H- x 1 
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Let x = a + ib, a, b real and presume b ^ 0 



We will show that this leads to a contradiction 



tan (a + ib) = (tana + i tanhb)/ ( 1 - i tana tanhb). 

Denoting the left side of (1) by L, then 



L = x tan x = (a + 
which may be written as 



ib j/ to* a + W i + j^fc lr 'ata~'£b \ 

\ I- \ / + {Tan a-la*fo) 



x tan x 



(ctA^ZCL- t> A^i-£zb)+i -tciA^-ixh ) 

Di 



( 2 ) 



(3) 



where 

= cos2a + cosh2b. 

Next considering the right hand side of equation (1) and denoting it 
by R we obtain 

^ _ C [(« -6-a z +h>-)-iiab>} (*-a l +b l ) + iLat> 

_ X ab (v'-a l -t-b l ) + 2.1L CL b 

, (4) 

(w-f- b*_ & l ) + 4# - & (>? V )- 2 l & 

— 

where 




a l J + 4a l b l 




(5) 



At this point it is well to observe that for b t 0, coshb ^ 1 

While for any real a, cosa > - 1 
and thus we conclude that 

D, > o (6' 
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Also 



t>i- 



(m-q + b 2 )-zi#b 



)> O H >o J b > , 



(7) 



2 7 

for (n - a +b ) and (2ab) cannot both vanish. 

Now if we equate the real and imaginary parts of (3) and (4), 
we obtain 



at - b/u^-dzb (*+b z -a 2 )-tqa l b L - £ (* + b z _ a) 



D, 



( 8 ) 



by equating the real parts, and 

/Wi 1 A Z b 

= _ z_£ 

t>. D i 



(9) 



by equating the imaginary parts. From (9) above we may write 



OL A*y\ 2. $ 

~dT 



q 1 2 b 

b Dj 



. l 

iSa. 

Dz 



(10) 



and substituting in (8) we obtain 

' b Dl t> z (U) 

Now considering the numerator on the right of (11). Call it N, and 
consider it as a function of S ; then 



— =- U + £-a) + l CL 2 = 3 a L -n-b 1 = f 

ds ^ 
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First, supposing that j 9 > 0, then the minimum value of N occurs 

for 6 = 0. 

Nmivt ” - a 2- ) + 4 a t > 0 

Next suppose that f <(. 0, then the minimum value of N occurs for 

& = n. 

+ k 2 - n ( * b 2 - ci 2 ,) + i n # 

1 (13) 

= n ( ol\ i> ) + b )> o 

Thus in any event /*/ ^ 0, and since for any real b ^ 0 we have 

> o (14) 

b 

thus the left side of (11) is less than zero whereas the right side is 
greater than zero; which is a contradiction. 

Thus we cannot have roots of the form x = a+ib with a and b real, 
and b ^ 0. Therefore we can only have real roots. 
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